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Efficient and stable algorithms for the calculation of spectral quantities and correlation functions 
are some of the key tools in computational condensed matter physics. In this article we review ba- 
sic properties and recent developments of Chebyshev expansion based algorithms and the Kernel 
Polynomial Method. Characterized by a resource consumption that scales linearly with the prob- 
lem dimension these methods enjoyed growing popularity over the last decade and found broad 
application not only in physics. Representative examples from the fields of disordered systems, 
strongly correlated electrons, electron-phonon interaction, and quantum spin systems we discuss 
in detail. In addition, we illustrate how the Kernel Polynomial Method is successfully embedded 
into other numerical techniques, such as Cluster Perturbation Theory or Monte Carlo simulation. 
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I. INTRODUCTION 

In most areas of physics the fundamental interactions 
and the equations of motion that govern the behavior of 
real systems on a microscopic scale are very well known, 
but when it comes to solving these equations they turn 
out to be exceedingly complicated. This holds, in par- 
ticular, if a large and realistic number of particles is in- 
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volved. Inventing and developing suitable approxima- 
tions and analytical tools has therefore always been a 
cornerstone of theoretical physics. Recently, however, 
research continued to focus on systems and materials, 
whose properties depend on the interplay of many dif- 
ferent degrees of freedom or on interactions that com- 
pete on similar energy scales. Analytical and approxi- 
mate methods quite often fail to describe the properties 
of such systems, so that the use of numerical methods re- 
mains the only way to proceed. On the other hand, the 
available computer power increased tremendously over 
the last decades, making direct simulations of the micro- 
scopic equations for reasonable system sizes or particle 
numbers more and more feasible. The success of such 
simulations, though, depends on the development and 
improvement of efficient algorithms. Corresponding re- 
search therefore plays an increasingly important role. 

On a microscopic level the behavior of most physical 
systems, like their thermodynamics or response to exter- 
nal probes, depends on the distribution of the eigenvalues 
and the properties of the eigenfunctions of a Hamilton op- 
erator or dynamical matrix. In numerical approaches the 
latter correspond to Hermitian matrices of finite dimen- 
sion D, which can become huge already for a moderate 
number of particles, lattice sites or grid points. The cal- 
culation of all eigenvalues and eigenvectors then easily 
turns into an intractable task, since for a D-dimensional 
matrix in general it requires memory of the order of D 2 , 
and the number of operations and the computation time 
scale as D 3 . Of course, this large resource consumption 
severely restricts the size of the systems that can be stud- 
ied by such a "naive" approach. For dense matrices the 
limit is currently of the order of D « 10 5 , and for sparse 
matrices the situation is only slightly better. 

Fortunately, alternatives are at hand: In the present 
article we review basic properties and recent develop- 
ments of numerical Chebyshev expansion and of the Ker- 
nel Polynomial Method (KPM). As the most time con- 
suming step these iterative approaches require only mul- 
tiplications of the considered matrix with a small set of 
vectors, and therefore allow for the calculation of spec- 
tral properties and dynamical correlation functions with 
a resource consumption that scales linearly with D for 
sparse matrices, or like D 2 otherwise. If the matrix is 
not stored but constructed on-the-fly dimensions of the 
order of D w 10 9 or more are accessible. 

The first step to achieve this favorable behavior is 
setting aside the requirement for a complete and exact 
knowledge of the spectrum. A natural approach, which 
has been considered from the early days of quantum me- 
chanics, is the characterization of the spectral density 
p{E) in terms of its moments fii = J p(E)E l dE. By 
iteration these moments can usually be calculated very 
efficiently, but practical implementations in the context 
of Gaussian quadrature showed that the reconstruction 
of p(E) from ordinary power mo ments is plagued by sub- 
stantial numerical instabilities ijGautschil Il968ft . These 
occur mainly because the powers E l put too much weight 



to the boundaries of the spectrum at the expense of poor 
precision for intermediate energies. The observation of 
this deficiency advanced the development of modified 
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mom ent approaches IjGautschil 1970; Sack and Donovan, 
11972ft . where E is replaced by (preferably orthogonal) 
polynomials of E. Wit h studies of the spectral den- 
sity of harmonic solids llBlumstein and Wheelerl Il973t 
IWheeler and Blumsteinlll972HWheeler et a/.lll974ft and 
of autocorrelation functions l)Wheelerlll974ft . which made 
use of Chebyshev polynomials of second kind, these ideas 
soon found their way into physics application. Later, 
similar Chebyshev expansion methods became popular 
also in quantum chemistry, where th e focus was on 
the time evolution of quantum states l| Chen an d GuoL 
[19991: iKoslofit Il988t iMandelshtam and Tavlorl Il997t 
iTal-Ezer an d Kosloff, 1984) and on Filter Diagonaliza- 
tion l)Neuhauserl ll99oTT ~ The modified moment ap- 
proach noticeably improved when kernel polynomials 
were introduced to damp the Gibbs oscillations, which 
for truncated polynomial series occ ur near disconti- 
nuitics of the expanded function (Silver and Rodcr, 
|l994 tlSilver et adll996tlWand.lT994llWang and Zuneerl 
1994). At this time also the name Kernel Polyno- 
mial Method was coined, and applications then in- 
cluded high- resolution spectral densities, static thermo- 
dynamic quantiti es as well as zero- t emperature dynam- 



ical correlations (Silv er and Rodeii 119941 IWanet fl99" 
IWang and Zungerl Il994ft . Only recently this range was 
extended to cover a lso dynamical correlation functions at 
finite-temperature (WeiBe, 2004), and below we present 
some new applications to complex-valued quantities, e.g. 
Green functions. Being such a general tool for studying 
large matrix problems, KPM can also be used as a core 
component of more involved numerical techniques. As re- 
cent examples we discuss Monte Carlo (MC) simulations 
and Cluster Perturbation Theory (CPT). 

In parallel to Chebyshev expansion techniques and 
to KPM also the Lancz o s Recursion Method was 
developed llAichhorn et al\ l2003t iBenoit et all 11992 ; 

iJaklic and Prelovsekl 



Havdock et all Il972l 



Lambin and Gaspardl 



1975 



1982) 



1994: 



whi ch is bas e d on a recur- 
sive Lanczos tridiagonalization l|Lanczosl Il950j) of the 
considered matrix and the expression of the spectral den- 
sity or of correlation functions in terms of continued frac- 
tions. The approach, in general, is applicable to the same 
problems as KPM and found wide application in solid 
state physics (|D ap;otto|. |l994l | Jak lic and Prelovsekl 12005 
lOrdeiorl Il998t iPantelidesl Il978ft . It suffers, however, 
from the shortcomings of the Lanczos algorithm, namely 
loss of orthogonality and spurious degeneracies if ex- 
tremal eigenstates start to converge. We will compare the 
two methods in Sec.lVland explain, why we prefer to use 
Lanczos for the calculation of extremal eigenstates and 
KPM for the calculation of spectral properties and corre- 
lation functions. In addition, we will comment on more 
specialize d itera t ive sc h emes, such as projection meth- 
ods llGoedeckerl Il999t iGoedecker and Colombo! Il994t 
llitaka and Ebisuzakil 12003ft and Maximum Entropy ap- 
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proaches llBandvopadhvav et all 120051: 1 Silver and Rodeit 
119971 ISkillineL Il988fl . Drawing more attention to KPM 
as a potent alternative to all these techniques is one of 
the purposes of the present work. 

The outline of the article is as follows: In Sec. |n] wc 
give a detailed introduction to Chebyshev expansion and 
the Kernel Polynomial Method, its mathematical back- 
ground, convergence properties and practical aspects of 
its implementation. In Sec. IIIII we apply KPM to a va- 
riety of problems from solid state physics. Thereby, we 
focus mainly on illustrating the types of quantities that 
can be calculated with KPM, rather than on the physics 
of the considered models. In Sec. II VI we show how KPM 
can be embedded into other numerical approaches that 
require knowledge of spectral properties or correlation 
functions, namely Monte Carlo simulation and Cluster 
Perturbation Theory. In Sec. El we shortly discuss al- 
ternatives to KPM and compare their performance and 
precision, before summarizing in Sec. IVII 



II. CHEBYSHEV EXPANSION AND THE KERNEL 
POLYNOMIAL METHOD (KPM) 

A. Basic features of Chebyshev expansion 

1. Chebyshev polynomials 

Let us first recall the basic properties of expansions in 
orthogonal polynomials and of Chebyshev expansion in 
particular. Given a positive weight function w(x) defined 
on the interval [a, b] we can introduce a scalar product 



(/Iff) = / w(x)f(x)g(x)dx 



(1) 



between two integrable functions /, g : [a, b] — » R. With 
respect to each such scalar product there exists a com- 
plete set of polynomials p n {x), which fulfil the orthogo- 
nality relations 



(2) 



where h n = l/(p n \p n ) denotes the inverse of the squared 
norm of p n (x) . These orthogonality relations allow for an 
easy expansion of a given function f{x) in terms of the 
p n (x), since the expansion coefficients are proportional 
to the scalar products of / and p n , 



fix) 



(x) with a n — (p n \f) h v 



(3) 



In general, all types of orthogonal polynomials can 
be used for such an expansion and for the Ker- 
nel Poly nomial approach we di scuss in this article 
(see e.g. ISilver and Roderl l)l994|) ). However, as we 
frequently observ e whenever we work with polyno- 
mial ^3xj}ans^ons_|Bo^|_^8^^ 

als ijAbramowitz and Stegunl Il970i iRivlinl Il990j) of first 



and second kind turn out to be the best choice for most 
applications, mainly due to the good convergence prop- 
erties of the corresponding series and to the close relation 
to Fourier transform llChene v Ul966tlLorentzUl966|) . The 
latter is also an important prerequisite for the derivation 
of optimal kernels (see Sec. I11.C|) . which are required for 
the regularization of finite-order expansions, and which 
so far have not been derived for other sets of orthogonal 
polynomials. 

Both sets of Chebyshev polynomials arc defined on the 
interval [a, 6] = [—1, 1], where the weight function w(x) — 
(71VI — x 2 )^ 1 yields the polynomials of first kind, T„, 
and the weight function w(x) = Try! — x 2 those of second 
kind, U n . Based on the scalar products 



, /|f/)i=/ M£K dx , 

7TV 1 — X Z 



(f\g)2= / TiVl -x 2 f(x)g(x)dx. 



the orthogonality relations thus read 



(T n |T m )i — 2 $n.m - 

2 

(U n \U m }2 — \ $n,m ■ 



(4) 



(5) 



(G) 
(7) 



By substituting x = cos(c/?) one can easily verify that they 
correspond to the orthogonality relations of trigonomet- 
ric functions, and that in terms of those the Chebyshev 
polynomials can be expressed in explicit form, 



T n (x) = cos(narccos(x)) , 

sin((n + 1) arccos(x)) 



U n {x) = 



sin(arccos(a;)) 



(8) 
(9) 



These expressions can then be used to prove the recursion 
relations, 



and 



T (x) = l, T_i(a;) = Ti(x) = x, 
T m+1 (x) = 2xT m (x) - T m _i(x) , 



U ix) = l, C/_!(x)=0, 
U m+ i(x) = 2xU m (x) - U m -i(x) , 



(10) 



(11) 



which illustrate that Eqs. © and © indeed describe 
polynomials, and which, moreover, are an integral part 
of the iterative numerical scheme we develop later on. 
Two other useful relations are 

2T m (x)T n (x) = T m+n (x) + T m _ n (x) , (12) 

2 (x 2 - 1) U m -i(x)U n -x(x) = T m+n (x) - T m ^ n (x) . 

(13) 

When calculating Green funct ions we also need Hilbert 
transforms of the polynomials l)Abramowitz and Stegunl 
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Il97d> . 



with moments 



F , T n (y)dy Un _ i(xh (M) 

[y - z)v 1 - v 



(f\<f>r, 



f(x)T n (x)dx . 



(23) 



p[ V^ U »-M* L = ^T n (x), (15) 

-l 

where V denotes the principal value. Chebyshev poly- 
nomials have many more interesting properties, for a de- 
tail ed discussion we refer the reader to text books such 
as [|Rivlinl.ll99oV 



2. Modified moments 

As sketched above, the standard way of expanding a 
function / : [—1,1] — ► R in terms of Chebyshev polyno- 
mials of first kind is given by 

f^) = T,WlTV T ^ x )= a o + 2j2^T n ( X ) (16) 



n=l 



with coefficients 



n„ = -/ir„ :i = / l^M dx . 

7TV 1 — X* 



(17) 



However, the calculation of these coefficients requires in- 
tegrations over the weight function w{x), which in prac- 
tical applications to matrix problems prohibits a simple 
iterative scheme. The solution to this problem follows 
from a slight rearrangement of the expansion, namely 



/(*) 



1 



7rVl — x 2 
with coefficients 



/i + 2^^„ T n (a 



fi n = / f(x)T n (x) dx . 



(18) 



(19) 



More formally this rearrangement of the Chebyshev series 
corresponds to using the second scalar product (.|.)2 and 
expanding in terms of the orthogonal functions 



<t>n{x) 



T n { x) 
ttvT" 



~2 ' 



which fulfil the orthogonality relations 

{4>n\4>m)2 = + 2'° 3n,m ■ 

The expansion in Eq. (|18fl is thus equivalent to 

fW = zrnor ^ 



(20) 
(21) 



ttVT 



Mo 



+ 2^ fi n T n (x) 



(22) 



The fj, n now have the form of modified moments that 
we announced in the introduction, and Eqs. (|18|l and l|19l) 
represent the elementary basis for the numerical method 
which we review in this article. In the remaining sections 
we will explain how to translate physical quantities into 
polynomial expansions of the form of Eq. (|18f) . how to 
calculate the moments /x ra in practice, and, most impor- 
tantly, how to regularize expansions of finite order. 

Naturally, the moments /j, n depend on the considered 
quantity f(x) and on the underlying model. We will 
specify these details when discussing particular applica- 
tions in Sec. II I II Nevertheless, there are features which 
are similar to all types of applications, and we start with 
presenting these general aspects in what follows. 



B. Calculation of moments 

1. General considerations 

A common feature of basically all Chebyshev expan- 
sions is the requirement for a rescaling of the underlying 
matrix or Hamiltonian H. As we described above, the 
Chebyshev polynomials of both first and second kind are 
defined on the real interval [—1, 1], whereas the quanti- 
ties we are interested in usually depend on the eigenval- 
ues {Ek} of the considered (finite-dimensional) matrix. 
To fit this spectrum into the interval [—1,1] we apply a 
simple linear transformation to the Hamiltonian and all 
energy scales, 



H = {H - b)/a. 
E = (E-b)/a, 



(24) 
(25) 



and denote all rescaled quantities with a tilde hereafter. 
Given the extremal eigenvalues of the Hamiltonian, E m [ n 
and -Emax, which can be ca lculated, e.g. with the Lanczos 
algorithm ijLanczosl ll95Cl(l . or for which bounds may be 
known analytically, the scaling factors a and b read 



a — (E nlax — E min )/(2 — e) , 
b = (E max + E m i n )/2 . 



(26) 
(27) 



The parameter e is a small cut-off introduced to avoid 
stability problems that arise if the spectrum includes or 
exceeds the boundaries of the interval [—1, 1]. It can be 
fixed, e.g. to e = 0.01, or adapted to the resolution of 
the calculation, which for an expansion of finite order N 
is proportional 1/N (see below). 

The next similarity of most Chebyshev expansions is 
the form of the moments, namely their dependence on 
the matrix or Hamiltonian H. In general, we find two 
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types of moments: Simple expectation values of Cheby- 
shev polynomials in H, 



Mr, 



'\T n (H)\a), 



(28) 



where \a) and \j3) are certain states of the system, or 
traces over such polynomials and a given operator A, 



M„ = TrL4T n (#)]. 



(29) 



Handling the first case is rather straightforward. Start- 
ing from the state \a) we can iteratively construct the 
states \a n ) = T n (H)\a) by using the recursion relations 
for the T n , Eq. JHJ, 



\a\) = H\a ) , 
\a n +i) = 2H\a n ) - \a n -i) ■ 

Scalar products with \(3) then directly yield 
fi„ = (P\a n ) . 



(30) 
(31) 
(32) 



(33) 



This iterative calculation of the moments, in particular 
the application of H to the state \a n ), represents the 
most time consuming part of the whole expansion ap- 
proach and determines its performance. If H is a sparse 
matrix of dimension D the matrix vector multiplication 
is an order 0(D) process and the calculation of N mo- 
ments therefore requires O(ND) operations and time. 
The memory consumption depends on the implementa- 
tion. For moderate problem dimension we can store the 
matrix and, in addition, need memory for two vectors 
of dimension D. For very large D the matrix certainly 
does not fit into the memory and has to be reconstructed 
on-the-fly in each iteration or retrieved from disc. The 
two vectors then determine the memory consumption of 
the calculation. Overall, the resource consumption of 
the moment iteration is similar or even slightly better 
than that of the Lanczos algorithm, which requires a few 
more vector operations (see our comparison in Sec. [VJ. 
In contrast to Lanczos, Chebyshev iteration is completely 
stable and can be carried out to arbitrary high order. 

The moment iteration can be simplified even further, 
if |/3) = | a). In this case the product relation l|12fl allows 
for the calculation of two moments from each new |a n ), 



\±2n = 2(a„|a„) - fi , 
H2n+i = 2(a n+ i|a„) - Hi 



(34) 
(35) 



which is equivalent to two moments per matrix vector 
multiplication. The numerical effort for N moments is 
thus reduced by a factor of two. In addition, like many 
other numerical approaches KPM benefits considerably 
from the use of symmetries that reduce the Hilbert space 
dimension. 



2. Stochastic evaluation of traces 

The second case where the moments depend on a trace 
over the whole Hilbert space, at first glance, looks far 
more complicated. Based on the previous considerations 
we would estimate the numerical effort to be proportional 
to D 2 , because the iteration needs to be repeated for all 
D states of a given basis. It turns out, however, that 
extremely good approximations of the moments can be 
obtained with a much sim pler approach: the stochas- 
tic evaluation of the trace llDrabold and Sankevt [l993: 
ISilver and Rodeil Il994t ISkillind. Il98^ ~ i.e., an estimate 
of fj, n based on the average over only a small number 
R <C D of randomly chosen states |r), 



Tr[AT n (H)] 



l^(r\AT n {H)\r) 



R 



(36) 



r=0 



The number of random states, R, does not scale with D. 
It can be kept constant or even reduced with increasing 
D. To understand this, let us consider the convergence 
properties of the above estimate. Given an arbitrary ba- 
sis {\i)} and a set of independent identically distributed 
random variables £h £ C, which in terms of the statisti- 
cal average ( • • ■ } fulfil 



«£ ri £r'i»=0, 

a random vector is defined through 

D-l 

\r) = Y,Ui)- 



(37) 
(38) 
(39) 



(40) 



We can now calculate the statistical expectation value of 

the trace estimate 0=4 J2^=o ( r \-^\ r ) f° r some Hermi- 
tian operator B with matrix elements Bij = (i\B\j), and 
indeed find, 



R-l 



R-l D-l 



r— r— i,j—0 

D-l 



i=0 



Of course, this only shows that we obtain the correct 
result on average. To assess the associated error we also 
need to study the fluctuation of O, which is characterized 
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by (SO) 2 = (e 2 )) - (e)) 2 . Evaluating 
«© 2 » = i^ E (r\B\r)(r'\B\r')} 

r,r'—0 

R-l D-l 
r,r'—0 . j' — 



r,r f —0 ij,i',j ,= 



R-l D-l 
r ,j'= 

(^ a + s(£W>*. 



D-l 



J 2 . 



3=0 



D-l 



E + E 



D-l 



(TrS) 2 + ^(Tr(B 2 ) + («|e„| 4 )) -2) £ S 



wc get for the fluctuation 



D-l 



(<5e) 2 = i(Tr(B 2 ) + ((|e„| 4 ))-2)^B 



3=0 



(42) 



(43) 



The trace of B 2 will usually be of order O(D), and the 
relative error of the trace estimate, 50/0, is thus of order 
0(1/ V RD). It is this favorable behavior, which ensures 
the convergence of the stochastic approach, and which 
was the basis for our initial statement that the number 
of random states R <C D can be kept small or even be 
reduced with the problem dimension D. 

Note also that the distribution of the elements of 
\r), p(£ri), has a slight influence on the precision of 
the estimate, since it determines the expectation value 
(l£™| 4 ) that enters Eq. fll3*)l. For an optimal distribution 
(I£h| 4 ) should be as close as possible to its lower bound 

(ICil 2 » = 1. and indeed > 

wc find this result if we fix the 
amplitude of the £ ri and allow only for a random phase 
4> G [0, 2tt], £ r i — e**. Moreover, if we were working in the 
eigenbasis of B this would cau s e (50 to vanish entirely, 
which led llitaka and Ebisuzakil l)2004Jl to conclude that 
random phase vectors are the optimal choice for stochas- 
tic trace estimates. However, all these considerations de- 
pend on the basis that we are working in, which in prac- 
tice will never be the eigenbasis of B (in particular, if B 
corresponds to something like AT n (H), as in Eq. JSfJJl). 
A random phase vector in one basis does not necessar- 
ily correspond to a random phase vector in another ba- 
sis, but the other basis may well lead to smaller value 
of Bjj , thus compensating for the larger value of 



(IC-il 4 )- Presumably, the most natural choice are Gaus- 
sian distributed which lead to (|£ ri | 4 )) = 2 and thus a 
basis-independent fluctuation (SO) 2 . To summarize this 
section, we think that the actual choice of the distribu- 
tion of £, ri is not of high practical significance, as long as 
Eqs. inn)-© are fulfilled for e C, or 



«eri»=Q, 



(44) 
(45) 



hold for £ ri E HL T ypically, within t his article we 
will consider Gaussian ijSilver and Ro deii ITofll ISklB 
Il988l) or uniformly distributed variables £ ri S R. 



C. Kernel polynomials and Gibbs oscillations 

1. Expansions of finite order & simple kernels 

In the preceding sections we introduced the basic ideas 
underlying the expansion of a function f(x) in an infinite 
series of Chebyshev polynomials, and gave a few hints for 
the numerical calculation of the expansion coefficients fj, n . 
As expected for a numerical approach, however, the total 
number of these moments will remain finite, and we thus 
arrive at a classical problem of approximation theory. 
Namely, we are looking for the best (uniform) approxi- 
mation to fix) by a polynomial of given maximal degree, 
which in our case is equivalent to finding the best approx- 
imation to f(x) given a finite number N of moments 
fi n . To our advantage, such problems have been stud- 
ied for at least 150 years and we can make use of results 
by many renowned mathematicians, such as Chebyshev, 
Weierstrass, Dirichlet, Fejer, Jackson, to name only a 
few. We will also introduce the concept of kernels, which 
facilitates the study of the convergence properties of the 
mapping f(x) — > /kpm(^) from the considered function 
f(x) to our approximation /kpm(i)- 

Experience shows that a simple truncation of an infi- 
nite series, 



i 



ttVI 



N-l 

! E 

71=1 



f^n T n (x) 



(46) 



leads to poor precision and fluctuations — also known 
as Gibbs oscillations — near points where the function 
f(x) is not continuously differentiable. The situation is 
even worse for discontinuities or singularities of f(x), as 
we illustrate below in Figure ^ A common procedure to 
damp these oscillations relies on an appropriate modifi- 
cation of the expansion coefficients, fj, n — > g n ^n, which 
depends on the order of the approximation TV, 

/KPMW = 2^ 9n JJ IA ; <Pn(x) 



71=0 



jV-1 



(47) 



ttVT 



go^o + 2 ^ 9nHn T n (x) 



n=l 
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In more abstract terms this truncation of the infinite se- 
ries to order N together with the corresponding modifi- 
cation of the coefficients is equivalent to the convolution 
of f(x) with a kernel of the form 



N-l 



K N (x,y) = g <l>o{x)(f>o(y) + 2 ^ g n <f>n(x)<f> n {y) , (48) 

71=1 

namely 



/kpm(z) = / Ti" V 1 -V 2 K N {x,y) f(y)dy 



(49) 



= (K N (x,y)\f(y))< 



The problem now translates into finding an optimal ker- 
nel K^{x,y), i.e., coefficients g n , where the notion of 
"optimal" partially depends on the considered applica- 
tion. 

The simplest kernel, which is usually attributed to 
Dirichlet, is obtained by setting g£ = 1 and evaluating 
the sum with the help of the Chris toffel-Darboux iden- 
tity ijAbramowitz and Stegurlll970|) . 

N-l 

Kn( x ,v) = <po(x)(/>o{y) + 2 ^2 <f>n(x)<l>n{y) 

«=i (50) 
_ (j>N{x)<j>N-i{y) - <pN-i(x)(l) N (y) 
x-y 

Obviously, convolution of with an integrable function 
/ yields the above truncated series, Eq. 1)46(1. which for 
N — > oo converges to / within the integral norm defined 
by the scalar product Eq. ||/||2 = \J (/|/)2, i-e. we 
have 



||/-/kpm|| 2 -5^0. 



(51) 



This is, of course, not particularly restrictive and leads 
to the disadvantages we mentioned earlier. 



2. Fejer kernel 



A first improvement is due to lFeierl l)l904[) who showed 
that for continuous functions an approximation based on 
the kernel 

1 N 

K§(x,y) = -J2^(x,y), i.e., ^ = 1-^,(52) 



converges uniformly in any restricted interval [— 1 + e, 1 — 
e]. This means that now the absolute difference between 
the function / and the approximation /kpm goes to zero, 



II/- /kpmII^o = , max \f{x) ~ /kpm(z)I N ^°°> 0. 

— l+e<ic<l — e 

(53) 



Owing to the denominator in the expansion (|46() con- 
vergence is not uniform in the vicinity of the endpoints 
x = ±1 , which we accounted for by the choice of a small 
e in the rescaling of the Hamiltonian H — > H . 

The more favorable uniform convergence is obtained 
under very general conditions. Specifically, it suffices to 
demand that: 

1. The kernel is positive: K?j(x,y) > OVx.y e [—1,1]. 

2. The kernel is normalized, f_..K(x,y)dx — (f>o(y), 
which is equivalent to go = 1. 

3. The second coefficient gi approaches 1 as N — > 00. 



Then , as a corollary to Korovkin's theorem (jKorovkinl 
1959), an approximation based on Kis;{x 1 y) converges 
uniformly in the sense explicated for the Fejer kernel. 
The coefficients g n , n > 2 are restricted only through 
the positivity of the kernel, the latter one being equiv- 
alent to monotonicity of the mapping / — > /kpm, i-e. 



/ > /' 



/] 



KPM 



> A 



KPM- 



Note also that the condi- 



tions ^ and |21 are very useful for practical applications: 
The first ensures that approximations of positive quan- 
tities become positive, the second conserves the integral 
of the expanded function, 



1 1 
JfKPM(x)dx = J f(x)dx. 

-1 -1 



(54) 



Applying the kernel, for example, to a density of states 
thus yields an approximation which is strictly positive 
and normalized. 

For a proof of the above theorem we refer th e reader 
to the literature llChenevl. Il96(t lLorentzl Il96fift . Let us 
here only check that the Fejer kernel indeed fulfils the 
conditions ^ to |21 The last two are obvious by inspection 
of Eq. (|52|l . To prove the positivity we start from the 
positive 27r-periodic function 



N-l 



(55) 



with arbitrary a„ € R. Straight-forward calculation then 
shows 



N-l 



N-l 



N-l N-lN-l-n 

= X] a " + 2 a vau+n cos rvp . 



v=0 



Hence, with 



n=l v=0 



N-l- 



(56) 



(57) 



v=0 
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the function 



and a„, respectively, note that 



N-l 



P(<P) = ,9o + 2 ^2 9n cosnip 



(58) 



is positive and periodic in (p. However, if p(ip) is positive, 
then the expression |[p(arccosx + arccosy)+p(arccosx — 
arccosy)] is positive V x, y € [— 1, 1]. Using Eq. (JSJ and 
cos a cos /3 = 2 [cos(a + (3) + cos(a — /?)], we immediately 
observe that the general kernel Kn(x, y) from Eq. is 
positive Va;,y S [—1,1], if the coefficients g n depend on 
arbitrary coefficients a v G M via Eq. Setting a„ = 

1/y/N yields the Fejer kernel K^(x, y), thus immediately 
proving its positivity. 

In terms of its analytical properties and of the conver- 
gence in the limit N — > oo the Fejer kernel is a major im- 
provement over the Dirichlet kernel. However, as yet we 
did not quantify the actual error of an order- N approxi- 
mation: For continuous functions an appropriate scale is 
given by the modulus of continuity, 



w f {A) 



max \f(x) 

v-y\<A 



f(y)\, 



in terms of which the Fejer approximation fulfils 



||/-/kpm|U~w/(1/VA0 



(59) 



(60) 



For sufficiently smooth functions this is equivalent to an 
error of order 0{l/yN). The latter is also an estimate for 
the resolution or broadening that we will observe when 
expanding less regular functions containing discontinu- 
ities or singularities, like the examples in Figure ^ 



3. Jackson kernel 

With the coefficients g£ of the Fejer kernel we have not 
fully exhausted the freedom offered by the coefficients a v 
and Eq. (|57|l . We can hope to further improve the kernel 
by optimizing the a„ in some sense, which will lead us to 
recover old results bv lJacksonl ljl91ll Il912|) . 

In particular, let us tighten the third of the previously 
defined conditions for uniform convergence by demanding 
that the kernel has optimal resolution in the sense that 



Q:= 



[x - y) 2 K N (x,y) dxdy 



(61) 



is minimal. Since Kn(x, y) will be peaked at x — y, Q is 
basically the squared width of this peak. For sufficiently 
smooth functions this more stringent condition will min- 
imize the error | \f — /kpm| |oo, and in all other cases lead 
to optimal resolution and smallest broadening of "sharp" 
features. 

To express the variance Q of the kernel in terms of g n 



( x — y) 2 = (Ti(x) — Ti(y)) 2 

= !(T 2 (a;) + T (x))T (y) - 2T 1 {x)T 1 {y) 

+ iT (x)(T 2 (y)+T (y)). (62) 

Using the orthogonality of the Chebyshev polynomials 
and inserting Eqs. I|48(l and l|62|l into (|f> 1 ft - we can thus 
rephrase the condition of optimal resolution as 



.9o — .9i = minimal w.r.t. a v . 



(63) 



Hence, compared to the previous section, where we 
merely required go — 1 and g± — > 1 for N — > oo, our 
new condition tries to optimize the rate at which g\ ap- 
proaches unity. 

Minimizing Q = g Q — gi under the constraint C = 
go — 1 = yields the condition 



dQ _ x dC 
da v da y 



(64) 



where A is a Lagrange multiplier. Using Eq. I|57|l and 
setting a_i = = we arrive at 



2a„ — a„_i — a u -\-\ — Aa„ 



(65) 



which the alert reader recognizes as the eigenvalue prob- 
lem of a harmonic chain with fixed boundary conditions. 
Its solution is given by 



a v = a sin 



A = 1 — cos ■ 



irk(v + 1) 
N+l 
n k 



(66) 



N + l 



where v = 0, . . . , (N - 1) and k = 1, 2, . . . , N. Given 
a v and the abbreviation q = nk/(N + 1) we can easily 
calculate the g n : 



2V-1- 



N- 



a v a v+n = a 2 sin qv sin q(v + n) 



y=0 
2 N—n 



= — [cos qn — cos q{2v + rij\ 



(67) 



a 

T 



JV- 



(N - n) cos qn - Re ^ e ifl< - 2u+ ^ 



= — [(N — n + l) cos qn + sin qn cot q] . 

The normalization go = 1 is ensured through a 2 — 
2/(N + 1), and with g x — cosq we can directly read off 
the optimal value for 



9o - ,9i = 1 - cos 



■nk 
N+l 



(68) 
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which is obtained for k = 1 

Qnun = 1 - COS - 



N 



7T 1 / 7T \ ' 



(69) 



The latter result shows that for large N the resolution 
•y/C? of the new kernel is proportional to 1/N. Clearly, 
this is an improvement over the Fejer kernel K^(x,y) 
which gives only ^[Q =1/ >/~N. 

With the above calc ulation we reproduced results 
bv lJacksonl (llQllL Il912ft . who showed that with a sim- 
ilar kernel a continuous function / can be approximated 
by a polynomial of degree N — 1 such that 



11/ -/kpm| \oo~W f (l/N). 



(70) 



which we may interpret as an error of the order of 
0(1/N). Hereafter we are thus referring to the new op- 
timal kernel as the Jackson kernel Kpj(x,y), with 

j (iV - n + 1) cos ^ +8^00^ 

9n N + 1 ■ K"-) 

Before proceeding with other kernels let us add a few 
more details on the resolution of the Jackson kernel: The 
quantity VQmin obtained in Eq. (j69|) is mainly a measure 
for the spread of the kernel K^(x,y) in the x-y-plane. 
However, for practical calculations, which may also in- 
volve singular functions, it is often reasonable to ask for 
the broadening of a <5-function under convolution with 
the kernel, 

*kpm(i - a) = (K N (x, y)\5(y ~ a)) 2 

N-l 

= go<j>o(x)T (a) + 2 ^ 9n4>n{x)T n {a) . (72) 
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Dirichlel kernel 
Jackson kernel 
Gaussian (a=ji/64) 
Lorentz kernel (X=4) 
Lorentzian (e=A/64) 





FIG. 1 (Color in online edition) Order N = 64 expansions of 
5(x) (left) and a step function (right) based on different ker- 
nels. Whereas the truncated series (Dirichlet kernel) strongly 
oscillate, the Jackson results smoothly converge to the ex- 
panded functions. The Lorentz kernel leads to relatively poor 
convergence at the boundaries x = ±1, but otherwise yields 
perfect Lorentz-broadened approximations. 



Using the Jackson kernel, an order TV expansion of a 5- 
function at x = thus results in a broadened peak of 
width a = whereas close to the boundaries, a = ±1, 



we hnd a 



N ' 



AT3/2 



It turns out that this peak is a good 



approximation to a Gaussian, 



V2t 



: exp i 



which we illustrate in Figure [I] 



4. Lorentz kernel 



It can be characterized by the variance a 1 = (x 2 )) — (x\ , 
where we use x = Ti(x) and x 2 — [T 2 (x) + T (x)]/2 to 
hnd 



{x} = J x5 K pm(x - a)dx = giTi(a) , (73) 

-l 
l 

g T (a) +g 2 T 2 (a) 



x 2 <5kpm(^ — a) dx 



(74) 



Hence, for K'^(x : y) the squared width of <5kpm(^ — a) is 
given by 



N -a 2 (N - 1) 
2(JV+1) 



I — cos ■ 



2tt 
7V + 1 



N 



3a 2 - 2 



N 



(75) 



The Jackson kernel derived in the preceding sections 
is the best choice for most of the applications we discuss 
below. In some situations, however, special analytical 
properties of the expanded functions become important, 
which only other kernels can account for. The Green 
functions that appear in the Cluster Perturbation The- 

Considering the imagi- 



ory, Sec. IIV.BI are an example, 
nary part of the Plemelj-Dirac formula which frequently 
occurs in connexion with Green functions, 



lim 



f 



+o x 



V 



i tt8(x) 



(77) 



the (5-function on the right hand side is approached in 
terms of a Lorentz curve, 



S(x) 



1 1 
- — lim Im 

7T e-+0 X + 



lim 



i e 



n(x 2 + e 2 ) 



(78) 



which has a different and broader shape compared to 
the approximations of S(x) we get with the Jackson 
kernel. There are attempts to approximate Lorentzian 
like behavior in the f ramework of filter diagonaliza- 
tion ijViiav et all |2004) , but these solutions do not lead 
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to a positive kernel. Note that positivity of the kernel 
is essential to guarantee basic properties of Green func- 
tions, e.g. that poles are located in the lower (upper) half 
complex plane for a retarded (advanced) Green function. 
Since we know that the Fourier transform of a Lorentz 
peak is given by exp(— e|fc|), we can try to construct an 
appropriate positive kernel assuming a„ = e~ Xl/ / N in 
Eq. I)57|). and indeed, after normalization, go = 1, this 
yields what we call the Lorentz kernel K^(x, y) hereafter, 



do much better, however, remembering the definition of 
the Chebyshev polynomials T n , Eq. JSJ, and the close 
relation between KPM and Fourier expansion: First, we 
may introduce the short-hand notation 

An = M«3n (81) 

for the kernel improved moments. Second and more im- 
portant, we make a special choice for our data points, 



L _ sinh[A(l - n/N)] 

Qn 



sinh(A) 



(79) 



The variable A is a free parameter of the kernel which 
as a compromise between good resolution and sufficient 
damping of the Gibbs oscillations we empirically choose 
to be of the order of 3 ... 5. It is related to the e- 
parameter of the Lorentz curve, i.e. to its resolution, via 
e = X/N. Note also, that in the limit A — > we recover 
the Fejer kernel K^(x, y) with = 1— n/N, suggesting 
that both kernels share many of their properties. 

In Figure ^ we compare truncated Chebyshev ex- 
pansions — equivalent to using the Dirichlet kernel - 
to the approximations obtained with the Jackson and 
Lorentz kernels, which we will later use almost exclu- 
sively. Clearly, both kernels yield much better approx- 
imations to the expanded functions and, in particular, 
the oscillations have disappeared almost completely. The 
comparison with a Gaussian or Lorentzian, respectively, 
illustrates the nature of the broadening of a ^-function 
under convolution with the kernels, which later on will fa- 
cilitate the interpretation of our numerical results. With 
Table [I] we conclude this section on kernels, and, for the 
sake of completeness, also list two other kernels that are 
occasionally used in the literature. Both have certain dis- 
advantages, in particular, they are not strictly positive. 



D. Implementational details and remarks 

1. Discrete cosine & Fourier transforms 

Having discussed the theory behind Chebyshev expan- 
sion, the calculation of moments, and the various kernel 
approximations, let us now come to the practical issues of 
the implementation of KPM, namely to the reconstruc- 
tion of the expanded function f{x) from its moments /i„. 
Knowing a finite number N of coefficients fj, n (see Sec. lIIII 
for examples and details), we usually want to reconstruct 
f(x) on a finite set of abscissas Xk- Naively we could sum 
up Eq. H47|) separately for each point, thereby making use 
of the recursion relations for T n , i.e., 



TTy/T 



5o Mo 



N-l 

71=1 



(80) 



For a set {xk} containing N points these summations 
would require of the order of NN operations. We can 



x k = cos v ' ' with k = 0, . . . , (N - 1) , (82) 

which coincides with the abscissas of Chebyshev numer- 
ical integration l|Abramowitz and Stegunl Il970l) . The 
number N of points in the set {x k } is not necessarily 
the same as the number of moments N. Usually we will 
consider N > N and a reasonable choice is, e.g. N = 2N. 
All values f(xk) can now be obtained through a discrete 
cosine transform, 



7fe = nyJl-X 2 k f(x k ) 



N-l 



MO 



2 ^ An cos ( 



n=l 



irn{k + 1/2) 
N 



(83) 



which allows for the use of divide-and-conquer type al- 
gorithms that require only N log N operations — a clear 
advantage over the above estimate NN. 

Routines for fast discrete cosine transform are im- 
plemented in many mathematical libraries or Fast 
Fourier Transform (FFT) packages, for instance, in 
FFTW ijFrigo and Johnsonl l2005albh that ships with 
most Linux distributions. If no direct implementation is 
at hand we may also use fast discrete Fourier transform. 
With 



'(2-<5„ )0 )A» exp(^) 



< n< N 
otherwise 



(84) 



and the standard definition of discrete Fourier transform, 



JV-l 



Afc = ^2 A„ exp 



n=0 



2nink 
N 



(85) 



after some reordering we find for an even number of data 
points 



j 2 j = Re(Aj) , 
72j+i = R«(A^_ 1 _ J .) . 



(86) 
(87) 



with j — 0, . . . , N/2 — 1. If we need only a discrete cosine 
transform this setup is not optimal, as it makes no use 
of the imaginary part which the complex FFT calculates. 
It turns out, however, that the "wasted" imaginary part 
is exactly what we need when we later calculate Green 
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Name 



fin 



Parameters positive? Remarks 



Jackson 
Lorentz 
Fcjcr 
Lanczos 

Wang and Zunger 
Dirichlet 



j^KN - n + 1) cos ^_ + sin ^_ cot ^] 
sinh[A(l -n/iV)]/sinh(A) 
l-n/N 

f sm(7rre/jV) \ A/ 
^ 7rn/JV I 

exp[-(ai)' 3 ] 



none yes best for most applications 

A € K yes best for Green functions 

none yes mainly of academic interest 

M 6 N no M = 3 closely matches the Jack- 

son kernel, but n ot strictly posi- 
tive iLanczosl ll966) 

a, /3 € R no found empirically, not op timal llWanel 
ITflflafeg and ZuweAlT^ 

none no least favorable choice 



TABLE I Summary of different integral kernels that can be used to improve the quality of an order TV Chebyshev series. The 
coefficients g n refer to Eq. 1471 or H48L respectively. 



(89) 



functions and other complex quantities, i.e., we can use 
the setup 



l2j = Aj , 

72J + 1 = Xfi-l-j > 

to evaluate Eq. (|14U[> . 

2. Integrals involving expanded functions 



We have already mentioned that our particular choice 
of Xk corresponds to the abscissas of numerical Cheby- 
shev integ ration. Hen c e, Ga uss-type numerical approx- 
imations ( Press et all . Il986|) to integrals of the form 
J , f{x)g(x)dx become simple sums, 



1 f{x)g{x)dx= j ^fmgjx) d:i 



VT^x* 



N-l 



N-l 



(90) 



N 



k=0 



k=0 



where 7^ denotes the raw output of the cosine or Fourier 
transforms defined in Eq. (|83(l . We can use this feature, 
for instance, to calculate partition functions, where f(x) 
corresponds to the expansion of the spectral density p(E) 
and g(x) to the Boltzmann or Fermi weight. 



E. Generalization to higher dimension 

1. Expansion of multivariate functions 

For the calculation of finite-temperature dynamical 
correlation functions we will later need expansions of 
functions of two variables. Let us therefore comment 
on the generalization of the previous considerations to d- 
dimensional space, which is easily obtained by extending 



the scalar product (.|.)2 to functions f,g : [—1, l] d — > '. 

1 1 d 
(f\gh = J ■ ■ J f(x)g(x) ( JJ - xf) dx! . . . dx 



-1 -1 



3=1 



(91) 

Here Xj denote the d components of the vector x. Natu- 
rally, this scalar product leads to the expansion 



m = E 



(fWh 

(4>n\4>7l) 
n=0 



(f>n(x) 



(92) 



where we introduced a vector notation for indices, n — 
{rix, . . . , nd}, and the following functions and coefficients 



Sri=0 Mra^n rij=l Irij {Xj) 



i>n(x) = JJ <Pnj(Xj) 
3=1 
= {f\<Pn)2 



(93) 



1 1 



ha 



/"' / II T nM)) dx X ... dx d , (94) 

-1 -1 ■?'=! 

1 d 2 



2. Kernels for multidimensional expansions 

As in the one-dimensional case, a simple truncation of 
the infinite series will lead to Gibbs oscillations and poor 
convergence. Fortunately, we can easily generalize our 
previous results for kernel approximations. In particular, 
we find that the extended kernel 



K N (x,y) = JJ K N (xj,yj) 
3=1 



(96) 
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maps an infinite series onto an truncated series, 
/kpm(z) = (K N (x,y)\f{y)} 2 



= T,nJ^ h " U a j=i 9n, T nA x j) (97) 



where we can take the g„ of any of the previously dis- 
cussed kernels. If we use the g J n of the Jackson kernel, 
K^(x,y) fulfils generalizations of our conditions for an 
optimal kernel, namely 

1. Kfj{x,y) is positive Vf,je [—1, l] d . 

2. K^(x,y) is normalized with 



ii 11 

• J /kpm(s) dxi . . . dx d = J ■■ J / (£) dx x ... dx d . 

l -l -l -l 

(98) 

3. Kfj(x,y) has optimal resolution in the sense that 
l l 

Q = J ■■ J(x - y) 2 K N (x, y) dx\ . . . dx d dy x ... dy d 



= d(g - .9i ) 
is minimal. 



(99) 



Note that for simplicity the order of the expansion, TV, 
was chosen to be the same for all spatial directions. Of 
course, we could also define more general kernels, 



K n ( f > V) = II Kn o ' % ) 

3=1 



(100) 



where the vector N denotes the orders of expansion for 
the different spatial directions. 



3. Reconstruction with cosine transforms 

Similar to the ID case we may consider the function 
/ : [— 1, 1] — ► R on a discrete grid Xt with 



x k,j = ^(Vk, ) , 

7r(fcj- + 1/2) 
N 

k, = 0, . . . , (N - 1) 



(101) 
(102) 
(103) 



Note again that we could define individual numbers of 

points for each spatial direction, i.e., a vector N with 
elements Nj instead of a single N. For all grid points 
x^ the function f(xt) is obtained through multidimen- 
sional discrete cosine transform, i.e., with coefficients 



«n = firthfi = Hftgnhft we find 

d 

Tfc = /(cos(^ fcl ), . . . ,cos(^ fc J) Y[ Trsm(ip k:i ) 

i=i 

N-l d 

= m kr n c ° s ( n j<Pk j ) (104) 

N-l N-l 
= 22 COS(ni99 fel )... ^ c °s( n dPk d )Kn ■ 
ni—0 ^d— 

The last line shows that the multidimensional discrete 
cosine transform is equivalent to a nesting of one- 
dimensional transforms in every coordinate. With fast 
implementations the computational effort is thus pro- 
portional to dN d ~ 1 N log N, which equals the expected 
value for N d data points, N d \ogN d . If we are not 
using libraries like FFTW, which provide ready-to-use 
multidimensional routines, we may also resort to one- 
dimensional cosine transform or the above translation 
into FFT to obtain high-performance implementations 
of general <i-dimensional transforms. 



III. APPLICATIONS OF KPM 

Having described the mathematical background and 
many details of the implementation of the Kernel Poly- 
nomial Method, we are now in the position to present 
practical applications of the approach. Already in the 
introduction we have mentioned that KPM can be used 
whenever we are interested in the spectral properties of 
large matrices or in correlation functions that can be ex- 
pressed through the eigenstates of such matrices. Appar- 
ently, this leads to a vast range of applications. In what 
follows, we try to cover all types of accessible quantities 
and for each give at least one example. We thereby focus 
on lattice models from solid state physics. 



A. Densities of states 

1. General considerations 

The first and basic application of Chebyshev expan- 
sion and KPM is the calculation of the spectral density of 
Hermitian matrices, which could correspond to the densi- 
ties of states of both interacting or non-interacting quan- 
tum models llSilver and Roderj 11994 ISilver et all 11996 " 



ISkillineL Il98a IWheelerlll974|) . To be specific, let us con- 
sider a Z)-dimensional matrix M with eigenvalues E/-, 
whose spectral density is defined as 



D-l 



P(E) 



5^ 



E k ) 



(105) 



fc=0 



As described earlier, the expansion of p{E) in terms of 
Chebyshev polynomials requires a rescaling of M — > M, 
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such that the spectrum of M — (M — b)/a lies within 
the interval [—1,1]. Given the eigenvalues of M the 
rescaled density p(E) reads 



D-l 
k=Q 



(106) 



and according to Eq. Ijl9(l the expansion coefficients be- 
come 



D-l 



Hn= f p(E) T n (E) dE=^J2 T «(^) 
= i E - T) Tr(T„(M)) 



(107) 



fe=0 



This is exactly the trace form that we introduced in 
Sec. III. Bl and we can immediately calculate the p n using 
the stochastic techniques described in Sec. I11.B.2I Know- 
ing the moments we can use the techniques of Sec. III.DI 
to reconstruct p(E) for the whole range [—1,1], and a 
final rescaling yields p(E). 



2. Non-interacting systems: Anderson model of disorder 

Applied to a generalized model of non-interacting 



fermions c. 



(t) 



D-l 
i,j=0 



(108) 



the matrix of interest M is formed by the coupling con- 
stants Mij. Knowing the spectrum of M, i.e. the single- 
particle density of states p{E) , all thermodynamic quan- 
tities of the model can be calculated. For example, the 
particle density is given by 



P{E) 



1 + e P( E -n) 
and the free energy per site reads 



dE 



(109) 



f = np- \ [ p{E) log(l + e - "^ - ^) dE , (110) 
P J 

where p is the chemical potential and (3=1 /T the inverse 
temperature. 

As the first physical example let us consider the An- 
derson model of non-interacting f ermions moving in a 
random potential l)Andersorlll958|) . 



(111) 



Here hopping occurs along nearest neighbor bonds 
(ij) on a simple cubic lattice and the local poten- 
tial ei is chosen randomly with uniform distribution 
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FIG. 2 (Color in online edition) Standard (dashed) and typi- 
cal density of states (solid line), p(E) and p typ (E) respectively, 
of the 3D Anderson model on a 50 3 site cluster with periodic 
boundary conditions. For p(E) we calculated = 2048 mo- 
ments with R — 10 start vectors and S = 240 realizations of 
disorder, for pt yp (E) these numbers are A^ = 8192, R = 32 
and S = 200. The lower right panel shows the phase diagram 
of the model we obtained from p typ (E)/ p(E) — ■> (mobility 
edge). 



in the interval [— W/2, W/2]. With increasing strength 
of disorder, W, the single-particle eigenstates of the 
model tend to become localized in the vicinity of 
a particular lattice site, which excludes these states 
from contributing to electronic transport. Disorder 
can therefore drive a transition from metallic behav- 
ior with delocalized f ermions to insulating behavior 
with localized fermions 1 Kramer and Ma c KinnonLll993t 
iLee and Ramakrishnanl 1985HThoulessl Il974|h The dis- 
order averaged density of states p(E) of the model can 
be obtained as described, but it contains no information 
about localization. The KPM method, however, allows 
also for the calculation of the local density of states, 



D-l 



M^) = ^Ei^ fc )i 2 ^-^)' ( 112 ) 

which is a measure for the contribution of a single lattice 
site (denoted by the basis state to the complete spec- 
trum. For delocalized states all sites contribute equally, 
whereas localized states reside on just a few sites, or, 
equivalently, a certain site contributes only to a few eigen- 
states. This property has a pronounced effect on the 
distribution of Pi(E), which at a fixed energy E char- 
acterizes the variation of pi over different realizations of 
disorder and sites i. For energies that correspond to lo- 
calized eigenstates the distribution is highly asymmet- 
ric and becomes singular in the thermodynamic limit, 
whereas in the delocalized case the distribution is regu- 
lar and centered near its expectation value p{E). There- 
fore a comparison of the geometric and the arithmetic 
average of p%{E) over a set of realizations of disorder 
and over lattic e sites reveals the position of the Ander- 
son transition ijDobrosavlievic and Kotliarl 119971 119981 
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ISchubert eitaZl . l2005albD . The expansion of Pi(E) is even 
simpler than the expansion of p(E), since the moments 
have the form of expectation values and do not involve a 
trace, 



} 1 

= \ Pi(E)T n {E)dE =pYl W\k)?T n (E k 

i k=0 



D-l 



(113) 



= £ £<«M)|fc)<fc|i) = i <i|T B (M)|t 



fe=0 



In Figure [21 we show the standard density of states p(E), 
which coincides with the arithmetic mean of Pi{E), in 
comparison to the typical density of states pt yp (E) , which 
is defined as the geometric mean of Pi(E), 



ftyp (£?)=ex P [(log(p i (£))>] 



(114) 



With increasing disorder, starting from the boundaries 
of the spectrum, p typ (E) is suppressed until it vanishes 
completely for W/t > 16.5, which is known as the critical 
strength of disord er where the states in the band center 
become localized (ISlevin and Ohtsukil [l999) . The calcu- 
lation yields the phase diagram shown in the lower right 
corner of Figure [21 which compares well to other numer- 
ical results. 

Since the method requires storage only for the sparse 
Hamiltonian matrix and for two vectors of the corre- 
sponding dimension, quite large systems can be stud- 
ied on standard desktop computers (of the order of 100 3 
sites). The recursion is stable for arbitrarily high expan- 
sion order. In the present case we calculated as many 
as 8192 moments to achieve maximum resolution in the 
local density of states. The standard density of states is 
usually far less demanding. 



3. Interacting systems: Double exchange 

Coming to interacting quantum systems, as a sec- 
ond example we study the evolution of the quantum 
double-exchange model ijAnderson and HaseeawaL Il955f) 
for large spin amplitude S, which in terms of spin-less 
fermions and Schwinger bosons a^J (a =1, |) is given 
by the Hamiltonian 



H = —- 



t 



2S+1 



a ia a ja c i c j 



(115) 



2S + c[ Cl . This 



with the local constraint Xo- a L a ier >w ~ " 
model describes itinerant electrons on a lattice whose 
spin is strongly coupled to local spins of amplitude S, 
so that the motion of the electrons mediates an effec- 
tive ferromagnetic interaction between these localized 
spin s. In the case of c olossal magneto-resistant mangan- 
ites l|Coev et adll999|) . for instance, cubic site symmetry 
leads to a crystal field splitting of the manganese d-shell, 



and three electrons in the resulting t2g-shell form the lo- 
cal spins. The remaining electrons occupy the e g -shell 
and can become itinerant upon dopin g, cau s ing t hese 
materials to show ferromagnetic order l)Zenerl Il95l|) . If 
the ferromagnetic (Hund's rule) coupling is large, at each 
site only the high-spin states are relevant and we can de- 
scribe the total on-site spin in terms of Schwinger bosons 
a^J l)Auerbachl . Il994j) . From the electrons only the 
charge degree of freedom remains, which is denoted by 
the spin-less fermions c-^ (see, e.g. ijWeifie et all l200lj) 
for more details). The full quantum model, Eq. (|115fl . 
is rather complicated for analytical or numerical studies, 
and we expect major simplification by treating the spin 
background classically (remember that S is quite large 
for the systems of interest). The limit of classical spins, 
S —>■ oo, is obtained by averaging Eq. (|1 15|) over spin 
coherent states, 

(cos(|) e^/ 2 a\ + sin(|) e" '*/ 2 a\) 2S 
\Q(S, 9, </>)) = ^ ^ T - I 2 ' y— |0) , 

(116) 

where 9 and </> are the classical polar angles and |0) the 
bosonic vacuum. The resulting non-interacting Hamilto- 
nian reads, 



H = - IZvls- + Hx - > 



(117) 



with the matrix element ijKogan and Auslenderl . fl9881 

t .. =i [cos|cos| e-^-tM 2 

+ sin | sin ^ e ^-«/ 2 ] , (118) 

i.e., spin-less fermions move in a background of random 
or ordered classical spins which affect their hopping am- 
plitude. 

To assess the quality of this classical approximation 
we considered four electrons moving on a ring of eight 
sites, and compared the densities of states obtained for 
a background of 5* = 3/2 quantum spins and a back- 
ground of classical spins. For the full quantum Hamil- 
tonian, Eq. I|115|) . the (canonical) density of states was 
calculated on the basis of 400 Chebyshev moments. To 
reduce the Hilbert space dimension and to save resources 
we made use of the SU (2) symmetry of the model: With 
the stochastic approach we calculated separate moments 



• '' for each S^-sector, 



= Tr^ [T n (H)} 



(119) 



and used the dimensions D s of the sectors to obtain the 
total normalized p, n from the average 



s ma3 

1 E 

Mn = ^Tr[T„(#)] = ^| 



,s* 



(120) 
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FIG. 3 (Color in online edition) Density of nonzero eigen- 
values of the quantum double-exchange model with S = 3/2 
(dashed line) and running average (red dot-dashed), calcu- 
lated for 4 electrons on a 8-site ring, compared to the classical 
result S — > oo (green solid). Expansion parameters: N — 400 
moments and R — 100 random vectors per S z sector. 



erator A reads 

(A) =^Tr(Ae-P H ) = ^£{k\A\k) e"^ , 

(121) 

z 4^ fl) 4^ e *' (i22) 

where 7? is the Hamiltonian of the system, Z the par- 
tition function, and Ek the energy of the eigenstate | fc) . 
Using the function 

a ( E ) = ^z2^ A ^ s ( E - E k) ( 123 ) 

fc=0 

and the (canonical) density of states p(E), we can express 
the thermal expectation value in terms of integrals over 
the Boltzmann weight, 



Note, that such a setup can be used whenever the model 
under consideration has certain symmetries. 

On the other hand, we solved the effective non- 
interacting model l|117ll and calculated the distributions 
of non-zero energies for a background of fully disordered 
classical spins. As Figure [21 illustrates, the spectrum of 
the quantum model with S = 3/2 closely matches that 
of the system with classical spins, providing good jus- 
tification, e.g. for studies of colossal magneto-resistive 
manganites that make use of a classical approximation 
for the spin background. Since for the finite cluster con- 
sidered the spectrum of the quantum model is discrete, 
at the present expansion order KPM starts to resolve 
distinct energy levels (dashed line). Therefore a running 
average (dot-dashed line) compares better to the classical 
spin-averaged data (bold line). 



(A) = - / a(E)e-' m dE, (124) 



Z= p(E) e~ pE dE. (125) 



Of course, similar relations hold also for non-interacting 
fermion systems, where the Boltzmann weig ht e~P E has 
to be replaced by the Fermi function f(E) — (1 + 
e /3(-E-A0)-r anc i |- ne si n gl e _electron wave functions play 
the role of | k) . 

Again, the particular form of a{E) suggests an expan- 
sion in Chebyshev polynomials, and after rescaling we 
find 

j D-l 

p„ = / a(E) T n (E) dE = - ^ (k\ A\k) T n {E k ) 

l x k=o (126) 

= ±-Tr(AT n (H)) , 



B. Static correlations at finite temperature 

Densities of states provide only the most basic infor- 
mation about a given quantum system, and much more 
details can usually be learned from the study of corre- 
lations and the response of the system to an external 
probe or perturbation. Starting with static correlation 
functions, let us now extend the application range of the 
expansion techniques to such more involved quantities. 

Given the eigenstates \k) of an interacting quantum 
system the thermodynamic expectation value of an op- 



which can be evaluated with the stochastic approach, 

Sec. EES 

For interacting systems at low temperature the expres- 
sion in Eq. (|124fl is a bit problematic, since the Boltz- 
mann factor puts most of the weight on the lower end 
of the spectrum and heavily amplifies small numerical 
errors in p(E) and a(E). We can avoid these problems 
by calculating the ground state and some of the lowest 
excitations exactly, using standard iterative diagonaliza- 
tion methods like Lanczos or Jacobi-Davidson. Then we 
split the expectation value of A and the partition func- 
tion Z into contributions from the exactly known states 
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FIG. 4 (Color in online edition) Nearest-neighbor S z -S z cor- 
relations of the XXZ model on a square lattice. Lines repre- 
sent the KPM results with separation of low-lying eigenstates 
(bold solid and bold dashed) and without (thin dashed), open 
symbols denote exact results from a complete diagonalization 
of a 4 x 4 system. 



and contributions from the rest of the spectrum 

c-i 



— OG 

(127) 



fe=0 



C-l 

k=0 J 



The functions 



D-l 



P> 



(E) = ±*£{k\A\k)S(E-E k ), 

k=C 

(E) = ^£5{E-E k ) 



(128) 



(129) 



(130) 



k=C 



describe the rest of the spectrum and can be expanded 
in Chebyshev polynomials easily. Based on the known 
states we can introduce the projection operator 

C-i 

P = l-^2\k)(k\, (131) 

fc=0 

and find for the expansion coefficients of a s {E) 

1 1 R_1 

(x n = -Tv(PAT n (H)) « — ^2(r\PAT n (S)P\r) , 



r=0 



(132) 



and similarly for those of p s (E) 

1 1 R ^ 

//„ = - Tv(PT n (H)) a — (r\PT n (H)P\r) . (133) 



Note, that in addition to the two vectors for the Cheby- 
shev recursion we now need memory also for the eigen- 
states \k). Otherwise the resource consumption is the 
same as in the standard scheme. 

We illustrate the accuracy of this approach in Figure^] 
considering the nearest-neighbor S z -S z correlations of 
the square-lattice spin- 1/2 XXZ model as an example, 



(134) 



i,<5 



RD 



r=0 



As a function of temperature and for an anisotropy 
— 1 < A < this model shows a quantum to clas- 
sical crossover in the sense that the correlations are 
anti-ferromagnetic at low temperature (quantum effect) 
and ferromagnetic at hig h temperature (as expected 
for the class i cal model) . llFabricius and McCovl Il999t 
iFehske et all 120001 ISchindelin et all l2000j) Comparing 
the KPM results with the exact correlations of a 4 x 4 
system, which were obtained from a complete diagonal- 
ization of the Hamiltonian, the improvement due to the 
separation of only a few low- lying eigenstates is obvious. 
Whereas for C — the data is more or less random below 
T w 1, the agreement with the exact data is perfect, if the 
ground state and one or two excitations are considered 
separately. The numerical effort required for these calcu- 
lations differs largely between complete diagonalization 
and the KPM method. For the former, 18 or 20 sites are 
practically the limit, whereas the latter can easily handle 
30 sites or more. 

Note that for non-interacting systems the above sepa- 
ration of the spectrum is not required, since for T — > 
the Fermi function converges to a simple step function 
without causing any numerical problems. 



C. Dynamical correlations at zero temperature 

1. General considerations 

Having discussed simple expectation values and static 
correlations, the calculation of time dependent quantities 
is the natural next step in the study of complex quantum 
models. This is motivated also by many experimental se- 
tups, which probe the response of a physical system to 
time dependent external perturbations. Examples are in- 
elastic scattering experiments or measurements of trans- 
port coefficients. In the framework of linear response 
theory and the Kubo formalism the system's response 
is expressed in terms of dynamical correlation functions, 
which can also be calculated efficiently with Chebyshev 
expansion and KPM. Technically though, we need to dis- 
tinguish between two different situations: For interacting 
many-particle systems at zero temperature only matrix 
elements between the ground state and excited states 
contribute to a dynamical correlation function, whereas 
for interacting systems at finite temperature or for non- 
interacting systems with a finite particle density transi- 
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tions between all eigenstates — many-particle or single- 
particle, respectively — contribute. We therefore split 
the discussion of dynamical correlations into two sections, 
starting here with interacting many-particle systems at 
T = 0. 

Given two operators A and B a general dynamical cor- 
relation function can be defined through 



<A;£)± = hm (0L4 



:B\0) 



uj + ie =F H 



k=0 



uj + ieTEk 



where Ek is the energy of the many-particle eigenstate 
of the Hamiltonian H, |0) its ground state, and e > 0. 

If we assume that the product (0|A|fc)(fe|B|0) is real 
the imaginary part 



D-l 



k=0 



has a similar structure as, e.g., the local density of states 
in Eq. (|112|) . and in fact, with pi(E) we already calculated 
a dynamical correlation function. Hence, after rescaling 
the Hamiltonian H — > H and all energies u> — > uj we can 
proceed as usual and expand lm(A;B)^ in Chebyshev 
polynomials, 



lm(A; B) 



/io + 2^^„T„(w) . (137) 

n=l 



Again, the moments are obtained from expectation values 



transform gives 



D-l 



MAB)t = J^(0\A\k)(k\B\0) V (-^ 



k=0 



1 v f Im(A; B} Q , dJ = _ 2 g ^ 



7r ._/ o; — uj 
-i 



(139) 



where we used Eq. (|14|) . The full correlation function 



vT 



n=l 



[^n-l(w) + 



i T n (u>) 



vT 



vT 



2 ^] Mr, 

n=l 



exp(— i n arccosw) 



(140) 

can thus be reconstructed from the same moments \x n 
that we derived for its imaginary part, Eq. i|138[) . In con- 
trast to the real quantities we considered so far, the re- 
construction merely requires complex Fourier transform, 
see Eqs. (|88|l and l|89|l . If only the imaginary or real part 
of (A; B)^ is needed, a cosine or sine transform, respec- 
tively, is sufficient. 

Note again, that the calculation of dynamical correla- 
tion functions for non-interacting electron systems is not 
possible with the scheme discussed in this section, not 
even at zero temperature. At finite band filling (finite 
chemical potential) the ground state consists of a sum 
over occupied single-electron states, and dynamical cor- 
relation functions thus involve a double summation over 
matrix elements between all single-particle eigenstates, 
weighted by the Fermi function. Clearly, this is more 
complicated than Eq. (|135fl . and we postpone the discus- 
sion of this case to Sec. llll.l51 where we describe methods 
for dynamical correlation functions at finite temperature 
and — for the case of non-interacting electrons — finite 
density. 



fi n = i J lm(A; B)± T n (uj) duj = (0\AT n (TH)B\0) , 

(138) 

and for A ^ B' we can follow the scheme outlined in 
Eqs. (JSUJl to For A = B^ the calculation simplifies 

to the one in Eqs. I|34|) and (|35J) . now with B|0) as the 
starting vector. 

In many cases, especially for the spectral functions and 
optical conductivities studied below, only the imaginary 
part of (A; B)^ is of interest, and the above setup is all 
we need. Sometimes however — e.g., within the Clus- 
ter Perturbation Theory discussed in Sec. IIV.BI — also 
the real part of a general correlation function (A; B)^ 
is required. Fortunately it can be calculated with al- 
most no additional effort: The analytical properties of 
(A; B)^ arising from causality imply that its real part is 
fully determined by the imaginary part. Indeed a Hilbert 



2. One-particle spectral function 

An important example of a dynamical correlation func- 
tion is the (retarded) Green function in momentum space, 

GAk,uj) = ( %a ;clX + (c^c,J Z , (141) 
and the associated spectral function 

A a (k,uj) = --lm G a (k,uj) 

7T (142) 

= A+(k,uj)+A-(k,uj), 

which characterizes the electron absorption or emission 
of an interacting system. For instance, A~ can be mea- 
sured experimentally in angle resolved photo-emission 
spectroscopy (ARPES). 



As the first application, let us consider the 
one-dime nsional Holstein model for spinlcss 
fermions l|Holsteml IT959albJ). 



i 

- guo X>1 + bi)m +"oJ2 b l b i > ( 143 ) 

which is one of the basic models for the study of electron- 
lattice interaction. A single band of electrons is approx- 
imated by spinless fermions cj , the density of which 
couples to the local lattice distortion described by dis- 
persionless phonons b\ . At low fermion density, with 
increasing electron phonon interaction the charges get 
dressed by a surrounding lattice distortion and form new, 
heavy quasi-particles known as polarons. Eventually, for 
strong coupling the width of the corresponding band is 
suppressed exponentially, leading to a process called self- 
trapping. For a half-filled band, i.e., 0.5 fermions per site, 
the model allows for the study of quantum effects at the 
transition from a metal to a band (or Peierls) insulator, 
marked by the opening of a gap at the Fermi wave vector 
and the development of a matching lattice distortion. 

Since the Hamiltonian (|143|l involves bosonic degrees 
of freedom, the Hilbert space of even a finite system has 
infinite dimension. In practice, nevertheless, the con- 
tribution of highly excited phonon states is found to be 
negligible at low temperature or for the ground-state, and 
the system is well approximated by a truncated phonon 
space with J2 t b l b i < M feauml et all Il998|) . In ad- 
dition, the translational symmetry of the model can be 
used to reduce the Hilbert space dimension, and, more- 
over, the symmetric phonon mode with momentum q = 
can be excluded from the numerics: Since it couples to 
the total number of electrons, which is a conserve d quan- 
tity, i ts contribution can b e handled analytically (Robin, 
19971 ISvkora et all l2005lh Below we present results for 
a cluster size of L = 8 or 10, where a cut-off M = 24 
or 15, respectively, leads to truncation errors < 10~ 6 
for the ground-state energy. Alternatively, for one or 
two fermionic particles and low temperatures an opti- 
mized va riational bas i s can be constructed for infinite 
systems i|Bonca et all ^999) , which would also be suit- 
able for our numerical approach. 

In Figure|3]we present KPM data for the spectral func- 
tion of the spinless-fermion Holstein model and assess its 
quality by comparing with results from Quantum Monte 
Carlo (QMC) and Dynamic al Density M a trix R cnormal- 
ization Group (DDMRG) jJeckelmanrJ, |2002|) calcula- 
tions. Starting with the case of a single electron on a 
ten-site ring, Figure |3] (a) illustrates the presence of a 
narrow polaron band at the Fermi level and of a broad 
range of incoherent contributions to the spectral func- 
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FIG. 5 (Color in online edition) One-particle spectral func- 
tion and its integral for the Holstein model (a) on a 10-site ring 
with one electron, e p = g 2 coo = 2.0i, loo = 0.4t, and (b) on a 
8-site ring, band filling n — 0.5, e p = g 2 tuo = 1.6t, wo = O.lt. 
For comparison, in (a) the blue das hed lines represent Quan- 
tum Monte Carlo data at (3t = 8 fH ohenadler et alj, 2005), 
and green stars ind icate the position o f the polaron band in 
the infinite system iBonca et al 1 119991) . In (b) the blue and 
green curves denote res ults of Dynamical DMRG for t he same 
lattice size and T = JJeckelmann and Fehskel l2006lh 



tion, which in the spinless case reads 



A-(k,io)=Y,\(l,N c -l\c k \0,N c )\ 2 
i 

x S[uj + (E^n^! - Eo^nJ] (144) 
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and 



A+{k,e>) = Y,\{l,N e + l\c{\0,N e )f 



x S[u- (E hNe+1 -Eo, N .)]. (145) 

Here \l, N e ) denotes the Ith eigenstate with N e elec- 
trons and energy Ei^- The photo-emission part A~ 
reflects the Poisson-like phonon distribution of the po- 
laron ground state, whereas A + has most of its weight in 
the vicinity of the original free electron band. In terms of 
the overall shape and the integrated weight, both KPM 
and QMC agree very well. QMC, however, is not able to 
resolve all the narrow features of the spectral function, 
and the polaron band is hardly observable. Nevertheless, 
QMC has the advantage that larger systems can be stud- 
ied, in particular at finite temperature. As a guide to the 
eye we also show the position of the polaron band in the 
in finite system, which was calculated with the approach 
of iBonca et all ill999|) . In Figure (b) we consider the 
case of a half-filled band and strong electron-phonon cou- 
pling, where the system is in an insulating phase with an 
excitation gap at the Fermi momentum k — ±7r/2. Be- 
low and above the gap the spectrum is characterized by 
broad multi-phonon absorption. Compared to DDMRG, 
again KPM offers the better resolution and unfolds all 
the discrete phonon sidebands. Concerning numerical 
performance DDMRG has the advantage of a small opti- 
mized Hilbert space, which can be handled with standard 
workstations. However, the basis optimization is rather 
time consuming and, in addition, each frequency value 
uj requires a new simulation. The KPM calculations, on 
the other hand, involved matrix dimensions between 10 8 
and 10 10 , and we therefore used high-performance com- 
puters such as Hitachi SR8000-F1 or IBM p690 for the 
moment calculation. For the reconstruction of the spec- 
tra, of course, a desktop computer is sufficient. 



3. Optical conductivity 

The next example of a dynamical correlation function 
is the optical conductivity. Here the imaginary and real 
parts of our general correlation functions (A; B) u change 
their roles due to an additional frequency integration. 
The so-called regular contribution to the real part of the 
optical conductivity is thus given by, 

o^( u ) = - V \(k\J\0)\ 2 6(uj - (E k - E Q )) , (146) 
u> z — ' 

E k >E 

where the operator 

J = -i<Z^(4 CT a i+1 , CT -H.c.) (147) 

describes the current. After rescaling the energy and 
shifting the frequency, uj = uj + E , the sum can be ex- 
panded as described earlier, now with J|0) as the initial 
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FIG. 6 (Color in online edition) (a) The optical conductiv- 
ity <T rog (aj) and its integral S Teg (u>) for the Holstein Hubbard 
model at half-filling with different ratios of the Coulomb inter- 
action U to the electron-lattice coupling e p = g 2 ioo, loo = O.lt, 
and g 2 = 7. Black dotted lines denote excitations of the pure 
Hubbard model, (b) The one-particle spectral function at the 
transition point, i.e., for the same parameters as in the middle 
panel of (a). The system size is L — 8. 



state for the Chebyshev recursion. Back-scaling and di- 
viding by uj then yields the final result. 

In Figure we apply this setup to the Holstein Hub- 
bard model, which is the generalization of the Holstein 
model to true, spin-carrying electrons that interact via 
a screened Coulomb interaction, modelled by a Hubbard 
/7-term, 

H=-t ^(4,a c i+l,<7 + H - C -) + U J2 n ^ Ui l 
i,a i 

-guJoY,( b i +b >i° +u} °Yl b i b i- ( 148 ) 

2,(7 i 

For a half-filled band, which now denotes a density of one 
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FIG. 7 Spin structure factor at T = calculated for the 
model (11491 which aims at describi ng the magnetic com pound 
(VO) 2 P2 7 . For more details see (IWeifie et aZ.U 19991) . 



electron per site, the electronic properties of the model 
are governed by a competition of two insulating phases: a 
Peierls (or band) insulator caused by the electron-lattice 
interaction and a Mott (or correlated) insulator caused 
by the electron-electron interaction. Within the opti- 
cal conductivity both phases are signalled by an exci- 
tation gap, which closes at the transition between the 
two phases. We illustrate this behavior in Figure (a), 
showing cr reg (<j;) at strong electron-phonon coupling and 
for increasing U. The data for the one-particle spectral 
function in FigureEl(b) proves that simultaneously to the 
optical gap also the ch arge gap vanishes at the q uantum 
phase transition point ijFehske et all 120041 120021. 



4. Spin structure factor 

Apart from electron systems, of course, the KPM ap- 
proach works also for other quantum problems such as 
pure spin systems. To describe the excitation spec- 
trum and the magnetic properties of the compound 
(VO)2P2C>7, some years ago we proposed the 2D spin 
Hamiltonian l|Weifle et n/1 llQQflfl 



H — J;,^J(1 + 5(— l) l )Sij ■ Si + ij + J a 2^ ' 
+ Jx y~^0^2i,j • Sai+ij+i + Sai+ij ■ sWj+i) , (149) 

where Sy denote spin- 1/2 operators on a square lattice. 
With this model we aimed at explaining the observation 



of two bran ches of low- lying trip let excitations by neutron 
scattering llGarrett et all Il997|) , which was inconsistent 
with the then prevailing picture of (VO)2P2C>7 being a 
spin-ladder or alternating chain compound. 

Studying the low-energy physics of the model l|149|) 
the KPM approach can be used to calculate the spin 
structure factor and the integrated spectral weight, 

S(q,u) = YsKkl^imfSiEk -Eo-u), (150) 

k 

N{q,w) = / du'S(q,uj'), (151) 
Jo 

where S z (q) — J2i j e lq ' Ti >i Sf j. Figure [7| shows these 
quantities for a 4 x 8 cluster with periodic boundary con- 
ditions. The dimension of the sector S z = 0, which con- 
tains the ground state, is rather moderate here being of 
the order of D w 4 • 10 7 only. The expansion clearly re- 
solves the lowest (massive) triplet excitations Ti , a num- 
ber of singlets and, in particular, a second triplet branch 
T%. The shaded region marks the two-particle contin- 
uum obtained by exciting two of the elementary triplets 
Ti, and illustrates that T2 is lower in energy. Since the 
system is finite in size, of course, the continuum appears 
only as a set of broad discrete peaks, the density of which 
increases with the system size. 



D. Dynamical correlations at finite temperature 

1. General considerations 

In the preceding section we mentioned briefly that for 
non-interacting electron systems or for interacting sys- 
tems at finite temperature the calculation of dynamical 
correlation functions is more involved, due to the required 
double summation over all matrix elements of the mea- 
sured operators. Chebyshev expansion, nevertheless, of- 
fers an efficient way for handling these problems. To 
be specific, let us derive all new ideas on the basis of 
the optical conductivity cr(u>), which will be our primary 
application below. Generalizations to other dynamical 
correlations can be derived without much effort. 

For an interacting system the extension of Eq. (|140H is 
given by 



E 

k.q 



\(k\J\q)\ 2 (e 



-PE k 



-PE q 



ZDuj 



S(U - UJ qk ) ■ 



(152) 

with uj q k — Eg — Ek- Compared to Eq. 1 146(1 a straight- 
forward expansion of the finite temperature conductiv- 
ity is spoiled by the presen ce of the Boltzmann w eight - 
ing factors. Some authors l)litaka and EbisuzakiL |2003|) 
try to handle this problem by expanding these factors 
in Chebyshev polynomials and performing a numerical 
time evolution subsequently, which, however, requires a 
new simulation for each temperature. A much simpler 
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(x,y) = ^^2\(k\J\q)\ 2 5{x-E k ) 5{y-E q ) (153) 



approach is based on the function 
1 

k,q 

which we may interpret as a matrix element density. Be- 
ing a function of two variables, j(x,y) can be expanded 
with two-dimensional KPM, 

j(x, y )= ^ — 2 m (154) 

where j{x,y) refers to the rescaled j{x,y), g n are the 
usual kernel damping factors (see Eq. I|71l0 . and h nm 
account for the correct normalization (see Eq. I|95[) ). The 
moments fi n m are obtained from 

l l 

/W = J J j(x,y)T n (x)T m (y)dxdy 
-l -l 



l 3 Y,MJ\q)\ 2 T n {E k ) T m {E q ) 

k.q 

3 Y.^\ T ^)J\q){q\T m {H)J\k) 



(155) 



k . 



= - Tr (T n (H) JT m (H) J) , 

and again the trace can be replaced by an average over a 
relatively small number R of random vectors |r). The 
numerical effort for an expansion of order n, m < N 
ranges between 2RDN and RDN 2 operations, depend- 
ing on whether memory is available for up to TV vectors 
of the Hilbert space dimension D or not. Given the op- 
erator density j(x, y) we find the optical conductivity by 
integrating over Boltzmann factors, 



oo 

^ rcg H = ^; J j(y + ",y)(< 



-0v . 



dy 
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-PE q 
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ZDuj 



■ S(UJ - LU qk ) 



(156) 

and, as above, we get the partition function Z from an 
integral over the density of states p(E). The latter can 
be expanded in parallel to j(x,y). Note that the cal- 
culation of the conductivity at different temperatures is 
based on the same operator density j(x,y), i.e., it needs 
to be expanded only once for all temperatures. 

Surprisingly, the basic steps of t his approach 
were suggested a l ready ten years ago ijWanet 1 1994 
IWang and Zungerl Il994|) . but — probably overlooking 
its potential — applied only to the zero-temperature re- 
sponse of non-interacting electrons. A reason for the poor 
appreciation of these old ideas may also lie in the use of 
non-optimal kernels, which did not ensure the positivity 




FIG. 8 (Color in online edition) The matrix element density 
j(x, y) for the 3D Anderson model with disorder W/t = 2 and 
12. 



of j(x,y) and reduced the numerical precision. Only re- 
cently, one of the authors generalized the Jackson kernel 
and obtain ed high resolu tion optical data for the Ander- 
son model (|WeifieL feof)^ . More results, in particular for 
interacting quantum systems at finite temperature, we 
present hereafter. 



2. Optical conductivity of the Anderson model 

Since the Anderson model describes non-interacting 
fcrmions, the eigenstates |fc) occurring in cr(a>) now de- 
note single-particle wave functions and the Boltzmann 
weight has to be replaced by the Fermi function, 

oo 

o**( w ) = I f j(y + u,y)(f(y)-f(y + u))dy 



^ \(k\J\q)\*(f(E k ) 

k.q 



f{E q )) 



5{w - uj qk ) 



(157) 

Clearly, from a computational point of view this expres- 
sion is of the same complexity for both, zero and finite 
temperature, and indeed, compared to Sec. IIII.CI we 
need the more advanced 2D KPM approach. 

Figure [S] shows the matrix element density j(x,y) cal- 
culated for the 3D Anderson model on a D = 50 3 site 
cluster. The expansion order is N — 64, and the mo- 
ment data was averaged over S = 10 disorder samples 
and R = 10 random start vectors each. Starting from a 
"shark fin" at weak disorder, with increasing W the den- 
sity j(x, y) spreads in the entire energy plane, simultane- 
ously developing a sharp dip along x — y. A comparison 
with Eq. i|157|) reveals that this dip is responsible for the 
decreas ing and final ly vanishing DC conductivity of the 
model llWeifiel 12004). In Figure we show the resulting 
optical conductivity at W/t — 12 for different chemical 
potentials \i and temperatures (3 = l/T. Note that all 
curves are derived from the same matrix element density 
j(x,y), which is now based onaD= 100 3 site cluster, 
expansion order N — 2048, an average over S = 440 
samples and only R — 1 random start vectors each. 
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FIG. 9 (Color in online edition) Optical conductivity of the 
3D Anderson model at disorder W — 12 and for different 
chemical potentials fi and temperatures /3 = 1/T. 



3. Optical conductivity of the Holstein model 

Having discussed dynamical correlations for non- 
interacting electrons, let us now come back to the case 
of interacting systems. The setup described so far works 
well for high temperatures, but as soon as T gets small we 
experience the same problems as with thermal expecta- 
tion values and static correlations. Again, the Boltzmann 
factors put most of the weight to the margins of the do- 
main of j(x,y), thus amplifying small numerical errors. 
To properly approach the limit T — > we therefore have 
to separate the ground state and a few lowest excita- 
tions from the rest of the spectrum in a fashion similar 
to the static correlations in Sec. IIII.BI Since we start 
from a 2D expansion, the correlation function (optical 
conductivity) now splits into three parts: a contribution 
from the transitions (or matrix elements) between the 
separated eigenstates, a sum of ID expansions for the 
transitions between the separated states and the rest of 
the spectrum (see Sec. 1111. CJl . and a 2D expansion for all 
transitions within the rest of the spectrum, 

C-l C-1D-1 D-l 

k,q=Q k=0 q=C k,q=C 



with 



Cfc, , 



\(k\J\q)\ 2 (e 



CT 2D (") 

(158) 



-0E h 



ZDlo 



(159) 



The expansions required for cr^f (o>) are carried out in 
analogy to Sec. IIII.C.3I but the resulting conductivities 
are weighted appropriately when all contributions are 
combined to cr reg (w). Using the projection operator de- 
fined in Eq. i"131|) . the corresponding moments read 
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FIG. 10 (Color in online edition) Left: Schematic setup for 
the calculation of finite-temperature dynamical correlations 
for interacting quantum systems, which requires a separation 
into parts handled by exact diagonalization (ED), ID Cheby- 
shev expansion and 2D Chebyshev expansion. Right: The 
lowest eigenvalues of the Holstein model on a six site chain 
for different electron-phonon coupling e p . The shaded region 
marks the lowest polaron band, which was handled separately 
when calculating the spectra in Figure [TT1 




FIG. 11 (Color in online edition) Finite temperature optical 
conductivity of a single electron coupled to the lattice via a 
Holstein type interaction. Different colors illustrate how, in 
particular, the low-temperature spectra ben efit from a sep- 
aratio n of C = 0, 1 or 6 low-energy states (ISchubert et all 
l2005d) . The phonon frequency is uo/t = 0.4. 



For cr%Y)(u)) we follow the scheme outlined in lIl~LD.il but 
use projected moments 



Mr, 



Tr{T n (H)PJT m (H)PJ)/D. 



(161) 



Mn = (k\JPT n (H)PJ\k) 



(160) 



In Figure IT0I we illustrate our setup schematically and 
show the lowest forty eigenvalues of the Holstein model, 
Eq. (|143(1 . with a band filling of one electron. Separating 
up to six states from the rest of the spectrum we obtain 
the finite-temperature optical conductivity of the system, 
Figure ITTI For high temperatures (T = t, see lower pan- 
els) the separation of low-energy states is not necessary, 
the conductivity curves for C = 0, 1 and 6 agree very 
well. For low temperatures (T = 0.1£, see upper panels), 
the separation is crucial. Without any separated states 
(C = 0) the conductivity has substantial numerical errors 
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and can even become negative, if large Boltzmann factors 
amplify infinitesimal numerical round-off errors of nega- 
tive sign. Splitting off the ground state (C = 1) or the 
entire (narrow) polaron band (C = 6 for the present six- 
site cluster), we obtain reliable, high-resolution spectra 
down to the lowest temperatures. From a physics point 
of view, at strong electron phonon coupling (right panels) 
the conductivity shows an interesting transfer of spectral 
weight from h igh to low frequencies, i f the temperature is 
increased (see ISchubert et al\ |2005c) for more details). 

With this discussion of optical conductivity as a finite 
temperature dynamical correlation function we conclude 
the section on direct applications of KPM. Of course, the 
described techniques can be used for the solution of many 
other interesting and numerically demanding problems, 
but an equally important field of applications emerges, 
when KPM is embedded into other numerical or analyt- 
ical techniques, which is the subject of the next section. 



IV. KPM AS A COMPONENT OF OTHER METHODS 
A. Monte Carlo simulations 

In condensed matter physics some of the most intensely 
studied materials are affected by a complex interplay of 
many degrees of freedom, and when deriving suitable 
approximate descriptions we frequently arrive at mod- 
els, where non-interacting fermions are coupled to classi- 
cal degrees of freedo m. Examples ar e colossal magneto- 
resistant ma neanites llDagottol [20031 or magnetic semi- 
conductors llSchliemanner^Il l20oTK . where the classi- 
cal variables correspond to localized spin degrees of free- 
dom. We already introduced such a model when we dis- 
cussed the limit 5* — > oo of the double-exchange model, 
Eq. (|117|l . The properties of these systems, e.g. a ferro- 
magnetic ordering as a function of temperature, can be 
studied by standard MC procedures. However, in con- 
trast to purely classical systems the energy of a given 
spin configuration, which enters the transition probabili- 
ties, cannot be calculated directly, but requires the solu- 
tion of the corresponding non-interacting fermion prob- 
lem. This is usually the most time consuming part, and 
an efficient MC algorithm should therefore evaluate the 
fermionic trace as fast and as seldom as possible. 

The first requirement can be matched by using KPM 
for calculating the density of states of the fermion sys- 
tem, which by integration over the Fermi function yields 
the energy of the underlying spin configuration. Com- 
bined with standard Metropolis single-spin updates this 
led to the first MC simulations of double - excha nge sys- 
tems (jMotome and FurukawaL Il999l 1200(1 l200l[) on rea- 
sonably large clusters (8 3 sites), which were later im- 
proved by replacing full traces by trace estimates and 
by increasi ng the efficiency of t he matrix vector multi- 
plicat ions ijAlvarez et all l2005t iFurukawa and Motomel 
12004ft . 

To fulfil the second requirement it would be advan- 
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FIG. 12 (Color in online edition) Magnetization as a func- 
tion of temperature for the classical double-exchange model 
at doping n — 0.5. We compare data obtained from 
the effective mo del iJ e ff (see text), from a hybrid Monte 
Carlo approach llAlonso et all 1200 if), the Truncated Poly- 
nomi al Expansion Method iMotome and Furukawal 120001 
I2001T) . and from a KPM based Cluster Monte Carlo tech- 
nique JWeifie et all 2005). L denotes the size of the under- 
lying three-dimensional cluster, i.e., D — L 3 is the dimension 
of the fermionic problem. 



tageous to replace the above single-spin updates by up- 
dates of the whole spin background. The first imple- 
mentation of such ideas was given in terms of an hybrid 
Monte Carlo algorithm ijAlonso et a^J . l200lf) . which com- 
bines an approximate time evolution of the spin system 
with a diagonalization of the fermionic problem by Leg- 
endre expansion, and requires a much smaller number of 
MC accept-reject steps. However, this approach has the 
drawback of involving a molecular dynamics type simu- 
lation of the classical degrees of freedom, which is a bit 
complicated and which may bias the system in the direc- 
tion of the assumed approximate dynamics. 

Focussing on the problem of classical double ex- 
change, Eq. 111711 . we the refore proposed a third ap- 
proach (|Weifie et all 120051. which combines the advan- 
tages of KPM w i th th e highly efficient Cluster M C al- 
gorithms llankel 119981: iKrautli 12004 IWolfl Il989h . In 
general, for a classical MC algorithm the transition prob- 
ability from state a to state b can be written as 



P{a 



A(a -> b)P(a -> b) , 



(162) 



where A{a — > b) is the probability of considering the move 
a — > 6, and P(a — > b) is the probability of accepting the 
move a — » b. Given the Boltzmann weights of the states 
a and b, W(a) and Wib), detailed balance requires that 



W(a)P{a -► b) = W(b)P(b -> a) , 



(163) 



which can be fulfilled with a generalized Metropolis al- 
gorithm 



P(a — > b) = min 1, 



W(b)A(b^a) 
W(a)A(a -> 6) 



(164) 
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In the standard MC approach for spin systems only a 
single randomly chosen spin is flipped. Hence, A(a — > 
b) = A(b — > a) and the probability P(a — ► 6) is usually 
much smaller than 1, since it depends on temperature 
via the weights W(a) and W(b). This disadvantage can 
be avoided by a clever construction of clusters of spins, 
which arc flipped simultaneously, such that the a priori 
probabilities A(a — > b) and A(b — ► a) soak up any differ- 
ence in the weights W(a) and W(b). We then ar rive at 
the famous rejection-free cluster MC algorithms (|Worft 
1989), which are characterized by P(a — * b) = 1. 

For the double-exchange model l|117|) we cannot expect 
to find an algorithm with P(a — > &) = 1, but even a 
method with P(a — * b) = 0.5 would be highly efficient. 
The amplitude of the hopping matrix element (|118fl is 
given by the cosine of half the relative angle between 
neighboring spins, or \Uj\ 2 = (1 + Si ■ Sj)/2. Averaging 
over the fermionic degrees of freedom, we thus arrive at 
an effective classical spin model 

H eS ^-Jes^^ + Si-Sj, (165) 

m 

where the particle density n approximately defines the 
coupling, J c ff n(l — n)/^/2. Similar to a classical 
Heisenberg model, the Hamiltonian H c g is a sum over 
contributions of single bonds, and we can therefore con- 
struct a cluster algorithm with P(a — > b) = 1. Surpris- 
ingly, the simulation of this pure spin model yields mag- 
netization data, which almost perfectly matches the re- 
sults for the full classical double-exchange model at dop- 
ing n = 0.5, see Figure 11*2*1 

For simulating the coupled spin fermion model l|117|l 
we suggested to apply the single cluster algorithm for 
H e g until approximately every spin in the system has 
been flipped once, thereby keeping track of all a priori 
probabilities A(a — > b) of subsequent cluster flips. Then 
for the new spin configuration the energy of the electron 
system is evaluated with the help of KPM. Note how- 
ever, that for a reliable discrimination of H c ff and the full 
fermionic model i|117|) the energy calculation needs to be 
very precise. For the moment calculation we therefore 
relied on complete trace summations instead of stochas- 
tic estimates. The KPM step is thus no longer linear 
in D, but still much faster than a full diagonalization 
of the bilinear fermionic model. Based on the resulting 
energy, the new spin configuration is accepted with the 
probability (|164f> . Figure []"""| shows the magnetization 
of the double-exchange model as a function of tempera- 
ture for n — 0.5. Except for small deviations near the 
critical temperature the data obtained with the new ap- 
proach co mpares well w i th the results of the hybrid MC 
approach ( Alo nso et all l200l|) . and due to the low nu- 
merical effort rather large systems can be studied. 

Of course, the combination of KPM and classical 
Monte Carlo not only works for spin systems. We may 
also think of models involving the coupling of electronic 
degrees of freedom to adiabatic lattice distortions or 



other classical variables ((Alvarez et all l2005[) , and as yet 
the potential of such combined approaches is certainly 
not fully exhausted. 

The next application, which makes use of KPM as a 
component of a more general numerical approach, brings 
us back to interacting quantum systems, in particular, 
correlated electron systems with strong local interactions. 

B. Cluster Perturbation Theory (CPT) 

1. General features of CPT 

Earlier in this review we have demonstrated the ad- 
vantages of the Chebyshev approach for the calculation 
of spectral functions, optical conductivities and struc- 
ture factors of complicated interacting quantum systems. 
However, owing to the finite size of the considered sys- 
tems, quantities like the spectral function A(k, lj) could 
only be calculated for a finite set of independent mo- 
menta k. The interpretation of this "discrete" data may 
sometimes be less convenient, e.g. the fc-integrated one- 
electron density p(w) — J dk d A(k,uj) does not show 
bands but only discrete poles which are grouped to band- 
like structures. Although this does not substantially bias 
the interpretation it is desirable to restore the transla- 
tional symmetry of the lattice and reintroduce an infinite 
momentum space. 

With the Cluster Perturbation T heory 
CCPT1 "Gros and Valentl Il994t ISenechal et all . l2000t 
ISenechal et 2Tl002) a straightforward way to perform 
this task approximatively has recently been devised. 
To describe it in a nutshell, let us consider a model of 
interacting fermions on a one-dimensional chain 

H = -tJ2(4+i,«^* + H.c.) + £ U t . (166) 

ia i 

Here f/j denotes a local interaction, e.g. [/, = Un^un 
for the Hubbard model. CPT starts by breaking up the 
infinite system into short finite chains of L sites each 
(clusters), which all are equivalent due to translational 
symmetry. From the Green function of a finite chain, 
Gij(aj) with i,j = 0,...,L—l, which is calculated ex- 
actly by a suitable numerical method, the Green function 
G(k,uj) of the infinite chain is obtained by reintroduc- 
ing the hopping between the segments. This inter-chain 
hopping is treated on the level of a random phase ap- 
proximation, which neglects correlations between differ- 
ent chains. The Green function G^ m (w) is then given 
through a Dyson equation 

i' ,j' ,m' 

(167) 

where V£ m = -t(5 n , m +i8ia8j,L-i + 5 n ,m-ih,L-i?>ja) de- 
scribes the inter-chain hopping and upper indices num- 
ber the different clusters. A partial Fourier transform 
of the inter-chain hopping, Vij(Q) = — ^e 1 ^ (^o^l-i + 
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FIG. 13 Spectral function for non-interacting tight-binding 
electrons. Based on the Lorentz kernel CPT exactly repro- 
duces the infinite system result (left), the Jackson kernel does 
not have the correct analytical properties, therefore CPT can- 
not close the finite size gap at k — it/2 (right). 



c l ® 5i t L-i6jo), gives the infinite-lattice Green function 
in a mixed representation 



G i3 (Q,uj) = 



1 - V{Q)G c (u;) 



(168) 



for a momentum vector Q of the super-lattice of finite 
chains and cluster indices Finally, from this mixed 
representation the infinite lattice Green function in mo- 
mentum space is recovered in the CPT approximation as 
a simple Fourier transform 



G(k,w) = j ^exp(i(i 



■j)k)Gij{Lk,uj). (169) 



i-,3 



The reader should be aware that restoring translational 
symmetry in the CPT sense is different from perform- 
ing the thermodynamic limit of the interacting system. 
The CPT may be understood as a kind of interpolation 
scheme from the discrete momentum space of a finite 
cluster to the continuous fc-values of the infinite lattice. 
The amount of information attainable from the solution 
of a finite cluster problem does however not increase. Es- 
pecially finite-size effects affecting the interaction prop- 
erties are by no means reduced, but still determined 
through the size of the underlying cluster. Nevertheless, 
CPT yields appealing presentations of the finite-cluster 
data, which can ease its interpretation. 

At present, all numerical studies within the CPT con- 
text use Lanczos recursion for the cluster diagonaliza- 
tion, thus suffering from the shortcomings we discussed 
earlier. As an alternative, we prefer to use the formalism 
introduced in Sec. IIII.CI which is much better suited for 
the calculation of spectral properties in a finite energy 
interval. 

On applying the CPT crucial attention has to be paid 
to the kernel used in the reconstruction of G^-(w). As 



it turns out, the Jackson kernel is an inadequate choice 
here, since already for the non-interacting tight-binding 
model it introduces spurious structures into the spectra. 
The failure can be attributed to the shape of the Jackson 
kernel: Being optimized for high resolution, a pole in the 
Green function will give a sharp peak with most of its 
weight concentrated at the center, and rapidly decaying 
tails. The reconstructed (cluster) Green function there- 
fore does not satisfy the correct analytical properties re- 
quired in the CPT step. To guarantee these properties, 
instead, we use the Lorentz kernel, which we constructed 
in Sec. III. C. 41 to mimic the effect of a finite imaginary 
part in the energy argument of a Green function. Us- 
ing this kernel for the reconstruction of G^AuS) the CPT 
works perfectly (cf. Figure ITS|l . 

To provide further examples we present results for two 
different interacting models where the cluster Green func- 
tion Gij(iv) has been calculated through a Chebyshev ex- 
pansion as in Eq. I|14U|) . Using G^Au) = G^(u;) (no mag- 
netic field), for a L-site chain L diagonal and L(L — l)/2 
off-diagonal elements of G^-(a->) have to be calculated. 
The latter can be reduced to Chebyshev iterations for 
the operators + , which allows application of the 
"doubling trick" (see the remark after Eq. (|138fl L How- 
ever, the numerical effort can be further reduced by a 
factor 1/L: If we keep the ground state |0) of the system 
we can calculate the moments (i% — (Q\ciT n {H)Cj\Q) for 
L elements i = 1, . ..,L of G^(uj) in a single Chebyshev 
iteration. To achieve a similar reduction within the Lanc- 
zos recursion we had to explicitly construct the eigen- 
states to the Lanczos eigenvalues. Then the factor 1/L 
is exceeded by at least ND additional operations for the 
construction of N eigenstates of a D-dimensional sparse 
matrix. Hence using KPM for the CPT cluster diagonal- 
ization the numerical effort can be reduced by a factor of 
1/L in comparison to the Lanczos recursion. 



2. CPT for the Hubbard model 

As a first example we consider the ID Hubbard model 
(Eq. I|148|) with q = uin = 0) , whic h is exactly solvable 
by Bethe ansatz l|Essler et adl2005l) and was also exten - 
sively studied with DDMRG ijJeckelmann et all l2000|) . 
It thus provides the opportunity to assess the precision 
of the KPM-based CPT. The top left panel of Figure HI 
shows the one-particle spectral function at half-filling, 
calculated on the basis of L = 16 site clusters and an 
expansion order of N = 2048. The matrix dimension is 
D sa 1.7 • 10 8 . Remember that the cluster Green function 
is calculated for a chain with open boundary conditions. 
The reduced symmetry compared to periodic boundary 
conditions results in a larger dimension of the Hilbcrt 
space that has to be dealt with numerically. In the top 
right panel the dots show the Bethe ansatz results for a 
L — 64 site chain, and the lines denote the L — > oo spinon 
and holon excitations each electron separates into (spin- 
charge separation). So far the Bethe ansatz does not 
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FIG. 14 (Color in online edition) Spectral function of the 
ID Hubbard model for half-filling and U = At. Top left: 
CPT result with cluster size L = 16 and expansion order 
N = 2048. For simila r data based on Lanczos recursion 
see lSenechal et al\ ||2000D . Top right: Within the exact Bethe 
ansatz solution each electron separates into the sum of inde- 
pendent spinon (red dashed) and holon (green) excitations. 
The dots mark the energies of a 64-site chain. Bottom: CPT 
data compared to selected DDMRG results for a system with 
L = 128 sites, open boundary conditions and a broadening of 
e = 0.0625£. Note that in DDMRG the momenta are approx- 
imate. 



allow for a direct calculation of the structure factor, the 
data thus represents only the position and density of the 
eigenstates, but is not weighted with the matrix elements 
of the operators cj^ . Although for an infinite system we 
would expect a continuous response, the CPT data shows 
some faint fine-structure. A comparison with the finite- 
size Bethe ansatz data suggests that these features are 
an artifact of the finite-cluster Greens function which the 
CPT spectral function is based on. The fine-structure is 
also evident in the lower panel of Figure 1141 where we 
compare with DDMRG data for a L = 128 site system. 
Otherwise the CPT nicely reproduces all expected fea- 
tures, like the excitation gap, the two pronounced spinon 
and holon branches, and the broad continuum. Note also, 
that CPT is applicable to all spatial dimensions, whereas 
DDMRG works well only for ID models. 
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FIG. 15 Spectral function A + (k, u) of a single electron in the 
Holstein model (corresponding to iV c = in Eq. 11451 '). For 
weak electron-phonon coupling the original band is still very 
pronounced (left), for intermediate-to-strong coupling many 
narrow polaron bands develop (right). The cluster size is 
L = 16 (l eft) or L = 6 (right) and the expansion order N = 
2048. SeelHohcnadl er~ef al\ J2003) for similar data based on 
Lanczos recursion. 

3. CPT for the Holstein model 

Our second example is the spectral function of a single 
electron in the Holstein model, i.e., Eq. (|148|) with U = 0. 
Here, as a function of the electron-phonon interaction, 
polaron formation sets in and the band width of the re- 
sulting quasi particles becomes extremely narrow at large 
coupling strength. Figure H31 illustrates this behavior for 
two values of the electron-phonon coupling e p = g 2 u>Q. 
For weak coupling the original one-electron band is still 
clearly visible (dot-dashed line), but the dispersion-less 
phonon (dashed line) cuts in approximately at an energy 
<x>o above the band minimum, causing the formation of 
a polaron band (solid line; calculated with the approach 
of lBonca et al\ l)l999jl ). an avoided-crossing like gap and 
a number of finite-size features. For strong coupling the 
spectral weight of the electron is distributed over many 
narrow polaron bands separated approximately by the 
bare phonon frequency luq. 

In all these cases, KPM works as a reliable high- 
resolution cluster solver, and using the concepts from 
Sec. llll.Dl we could also extend these calculations to finite 
temperature. Probably, CPT is not the only approximate 
technique that profits from the simplicity and stability of 
KPM, and the range of its applications can certainly be 
extended. 



V. KPM VERSUS OTHER NUMERICAL APPROACHES 

After we have given a very detailed description of the 
Kernel Polynomial Method and presented a wide range 
of applications, let us now classify the method in the 
context of numerical many-particle techniques and com- 
ment on a number of other numerical approaches that 
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are closely related to KPM. 



A. KPM and dedicated many-particle techniques 



In the previous sections we already compared KPM 
data and results of other numerical many-particle tech- 
niques. Nevertheless, it seems appropriate to add a 
few comments about the general concept of such cal- 
culations and the role KPM-like methods play in the 
field of many-particle physics and complex quantum sys- 
tems. The numerical study of interacting quantum many- 
particle systems is complicated by the huge Hilbcrt space 
dimensions involved, which usually grow exponentially 
with the number of particles or the system size. There 
are different strategies to cope with this: In Monte 
Carlo approaches only part of the Hilbert space is sam- 
pled stochastically, thereby trying to capture the es- 
sential physics with an appropriate weighting mecha- 
nism. On the other hand, variational metho ds, like 
DMRG llPeschel rt. all ISchollwockt Eool or the 

specialized approach of iBonca et all l)l999j) . aim at re- 
ducing the Hilbert space dimension in an intelligent way 
by discarding unimportant states, which, for instance, 
contribute only at high temperature. Compared to such 
methods KPM is much more basic: It is designed only for 
the fast and stable calculation of the spectral properties 
of a given matrix and of related correlations. Choosing a 
suitable Hilbert space or optimizing the basis is the mat- 
ter of the user or of external programs. It is thus a more 
general approach, which can be used directly or embed- 
ded into other methods, as we illustrated in the preced- 
ing section. Of course, this simplicity and general appli- 
cability come at a certain price: For interacting many- 
particle models the system sizes that can be studied by 
using KPM directly are usually much smaller, compared 
to DMRG and Monte Carlo. Note however, that both of 
the latter methods have limitations too: For many inter- 
esting models Monte Carlo methods are plagued by the 
infamous sign problem, which is not present in KPM. 
When it comes to the calculation of dynamical correla- 
tion functions Monte Carlo approaches rely on power mo- 
ments. The reconstruction of correlation functions from 
power moments is known to be an ill-conditioned prob- 
lem, in particular, if the moments are subject to sta- 
tistical noise. The resolution of Monte Carlo results is 
therefore much smaller compared to the data obtained 
with KPM. The DMRG method develops its full poten- 
tial only in one spatial dimension and for short ranged 
interactions. In addition, the calculation of dynamical 
correlations is lim ited to zero temperature, with only 
a few exceptions (jSirker and Klumperl . [2005) . None of 
these restrictions apply to KPM. 
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FIG. 16 (Color in online edition) Comparison of a KPM and a 
MEM approximation to a spectrum consisting of five isolated 
rj-peaks, and to a step function. The expansion order is N = 
512. Clearly, for the <5-peaks MEM yields a higher resolution, 
but for the step function the Gibbs oscillations return. 



B. Close relatives of KPM 

Having compared KPM to specialized many-particle 
methods, let us now discuss more direct competitors of 
KPM, i.e., methods that share the broad application 
range and some of its general concepts. 



1. Chebyshev expansion and Maximum Entropy Methods 

The first of these approaches, the combination of 
Chebyshev expansion and Maximum Entropy (MEM) , is 
basically an alternative procedure to transform moment 
data fi n into convergent approximations of the considered 
function f(x). To achieve this, instead of (or in addition 
to) applying kernel polynomials, an entropy 

S(f,fo)= J U(x)-f (x)-\og(f(x)/Mx)))dx (170) 

is maximized under the constraint that the moments of 
the estimated f{x) agree with the given data. The func- 
tion fo(x) describes our initial knowledge about /(x), 
and may in the worst case just be a constant. Being 
related to Maximum E ntropy approaches to th e clas- 
sical momen t problem l|Mead and Papanicolaou! Il984t 
iTurekl [ 1988), for the case of Chebyshev moments dif- 
ferent i mplementations of the meth o d have been sug- 
gested l| Bandyo padhvav et all 120051: I Silver and R5deit 
Since for a given set of N moments 
Hn the approximation to the function f{x) is usually not 
restricted to a polynomial of degree N — 1, compared 
to the KPM with Jackson kernel the Maximum Entropy 
approach usually yields estimates of higher resolution. 
However, this higher resolution results from adding a pri- 
ori assumptions and not from a true information gain (see 
also Figure The resource consumption of Maximum 
Entropy is generally much higher than the NlogN be- 
havior we found for KPM. In addition, the approach is 
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non-linear in the moments and can occasionally become 
unstable for large N. Note also that as yet Maximum En- 
tropy methods have been derived only for positive quan- 
tities, f(x) > 0, such as densities of states or strictly 
positive correlation functions. 

Maximum Entropy nevertheless, is a good alternative 
to KPM, if the calculation of the /i„ is particularly time 
consuming. Based on only a moderate number of mo- 
ments it yields very detailed approximations of /(x), and 
we obtained very goo d results for some c omputationally 
demanding problems (Bau ml et all 119981) . 



2. Lanczos recursion 

The Lanczos Recursion Method is certainly the 
most ca pable competit or of the Kernel Polynomial 
Method llDagottol Il994l) It is based on the Lanczos 
algorithrri ljLanczosl 11950^1 . a method which was ini- 
tially developed for the tridiagonalization of Hermitian 
matrices and later evolved to one of the most power- 
ful methods for th e calculation of extremal eigen states 
of sparse matrices l)Cullum and Willoughbvi Il985l) . Al- 
though ideas like the mapping of the classical moment 
problem to tridiagonal matri ces and continu ed fractions 
have been suggested earlier (iGordorl Il96l . the use of 
the Lanczos a lgorithm for the characteri zation of spec- 
tral densities IjHavdock et all . Il972t Il975|) was first pro- 
posed at about the same time as the Chebyshev expan- 
sion approaches, and in principle Lanczo s recursion is 
also a kind of modified mom e nt exp ansion ijBenoit et al\ . 
Il992t lLambin and Gasoardl Il982|) . Its generalization 
from spectral densities to zero temperature dynamical 
correlation fun ctions was first given in term s of contin- 
ued fractions ijGagliano and Balseirol Il987j) . and later 
also an approach based on the eigenstates of the tridiago- 
nal matri x was introduced an d termed Spectral Decoding 
Method ijZhong et all Il994l) . Thi s technique was then 
generalized to finite temperature IjJaklic and Prelovsekl 
1994 1200(1 . and, in addition, some variants of the 
approach for low temperature l)Aichhorn et_oj] 1 _|2003 ) 
and b ased on the micro-canonical ensemble l|Long"e^I , 
2003) have been proposed recently. 

To give an impression, in Table [n] we compare the 
setup for the calculation of a zero temperature dynamical 
correlation function within the Chebyshev and the Lanc- 
zos approach. The most time consuming step for both 
methods is the recursive construction of a set of vectors 
\4> n ) i which in terms of scalar products yield the moments 
fi n of the Chebyshev series or the elements a n , p n of the 
Lanczos tridiagonal matrix. In terms of the number of 
operations the Chebyshev recursion has a small advan- 
tage, but, of course, the application of the Hamiltonian 
as the dominant factor is the same for both methods. As 
a drawback, at high expansion order the Lanczos itera- 
tion tends to lose the orthogonality between the vectors 
1 4> n )i which it intends to establish by construction. When 
the Lanczos algorithm is applied to eigenvalue problems 



this loss of orthogonality usually signals the convergence 
of extremal eigenstates, and the algorithm then starts to 
generate artificial copies of the converged states. For the 
calculation of spectral densities or correlation functions 
this means that the information content of the a n and 
(3 n does no longer increase proportionally to the num- 
ber of iterations. Unfortunately, this deficiency can only 
be cured with more complex variants of the algorithm, 
which also increase the resource consumption. Cheby- 
shev expansion is free from such defects, as there is a 
priori no orthogonality between the \4>n)- 

The reconstruction of the considered function from its 
moments /i n or coefficients a„, /3„, respectively, is also 
faster and simpler within the KPM, as it makes use of 
Fast Fourier Transformation. In addition, the KPM is 
a linear transformation of the moments a property 
we used extensively above when averaging moment data 
instead of the corresponding functions. Continued frac- 
tions, in contrast, are non- linear in the coefficients a n , 
(3 n . A further advantage of KPM is our good under- 
standing of its convergence and resolution as a function 
of the expansion order N. For the Lanczos algorithm 
these issues have not been worked out with the same 
rigor. 

We therefore think that the Lanczos algorithm is an 
excellent tool for the calculation of extremal eigenstates 
of large sparse matrices, but for spectral densities and 
correlation functions the Kernel Polynomial Method is 
the better choice. Of course, the advantages of both al- 
gorithms can be combined, e.g. when the Chebyshev 
expansion starts from an exact eigenstate that was cal- 
culated with the Lanczos algorithm. 



3. Projection methods 

Projection methods were developed mainly in the con- 
text of electronic structure calculations or tight-binding 
molecular dynamics, which both require knowledge of 
the total energy of a non-int eracting electron system or 
of rel ated expectation values l)Goedeckeit Il999t lOrdeionl 
1998). The starting point of these methods is the den- 
sity matrix F = f(H), where f(E) again represents 
the Fermi function. Thermal expectation values, total 
energies and other quantities of interest are then ex- 
pressed in terms of traces over F and co rresponding op- 
erators ijGoedecker and Colombol Il994j) . For instance, 
the number of electrons and their energy are given by 
Aci = Tt(F) and E = Tr(FH), respectively. To obtain 
a numerical approach that is linear in the dimension D 
of H, F is expanded as a series of polynomials or other 
suitable functions in the Hamiltonian H, 

1 N ~ X 

J1 = 1+C ^-M) -E a ^( g )> ( m ) 

and the above traces are replaced by averages over ran- 
dom vectors \r). Chebyshev polynomials are a good basis 
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Chebyshev / KPM 



complexity Lanczos recursion 



complexity 



Initialization: 



H = (H-b)/a 
\<h)=A\0), \<h) = H\<h) 
jtto = (<£o |<M, Mi = (0i|0o) 



Initialization: 



/3 = v7oRtA|0)" 
\4>o)=A\0)/Po, |0-i) 



Recursion for 2A r moments fj, n : 

\4>n + l) = 2H\(p„) — \<j>n-l) 
fJ-2n + 2 = 2(0 n + i|0 n+ i) — fJ,Q 
M2n + 1 = 2(0 n +l|0„) — Hi 



O(ND) Recursion for N coefficients a n , Pn- 

\<j>) = H\(j)n) — f3n\(f>n-l), Ct„ = ((f>n\<t>') 



') = \<f>') -a n \4> n ), /?»+! = V^WO 
|0n+l> = \<t>")/Pn+l 



O(ND) 



very stable 



tends to lose orthogonality 



Reconstruction in three simple steps: 
Apply kernel: /i n = g n ^n 
Fourier transform: /Et n — * f(u>i) 
f[( Ui - b)/a] 



Rescale: 



0(M log M) Reconstruction via continued fraction 



f(z) = --Jm- 
n 



z — ao — 



fit 



Z — Oil — 



Z — Oil 



where z = Wj + i e 



O(NM) 



procedure is linear in /i n — > procedure is non-linear in a n , /3, 
well defined resolution oc 1/JV — > e is somewhat arbitrary 



TABLE II Comparison of Chebyshev expansion and Lanczos recursion for the calculation of a zero-temperature dynamical 
correlation function /(w) = ~^2 n \ (n\A\0) \ 2 S(u> — ui„). We assume N matrix vector multiplications with a D-dimensional sparse 
matrix H, and a reconstruction of f(u) at M points u)i. 



for such an expansion of F l)Goedecker and TeterL [199^) ■ 
and the corresponding approaches are thus closely related 
to the KPM setup we described in Sec. IIII.AI Note how- 
ever, that the expansion in Eq. 1|171|) has to be repeated 
whenever the temperature 1//3 or the chemical poten- 
tial \i is modified. This is particularly inconvenient, if n 
needs to be adjusted to fix the electron density of the sys- 
tem. To compensate for this drawback, at least partially, 
we can make use of the fact that in Eq. (|171fl the ex- 
panded function and its expansion coeffic ients are kn own 
in advance: Using implicit methods ijNiklassori |2003|) the 
order N approximation of F can be calculated with only 
(log TV) matrix vector operations involving the Hamil- 
tonian H . The total computation time for one expansion 
is thus proportional to DlogN, compared to DN if the 
sum in Eq. I|171|) is evaluated iteratively, e.g., on the basis 
of the recursion relation Eq. l|lUfl • 



Projection methods can also be used for the calcula- 
tion of dynamical correlation functions. In this case the 
expansion of the density matrix, which accounts for the 
thermodynamics, is supplemented by a numerical time 
evolution. Hence, a general correlation function is writ- 
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FIG. 17 (Color in online edition) The optical conductivity of 
the Anderson mode l. Eq. calculated with KPM and a 

projection method llitaka ll99Si) . The disorder is W ~ 15; 
temperature and chemical potential read T — and /i — 0. 
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ten as 

oo 

{A;B) U = Urn J e i(w+ie)t Tr(e [Ht Ae~ [Ht BF)dt, 
o 

(172) 

and the e ±lHt terms are handled b y standard 
meth ods, such as Crank- Nicolson lj£ress_et_al , 
1986), Suzuki- Trotter ijde Vries and De Raedi . 
1993|) , and, very efficiently, Che byshev expan- 
sion ijDobrovitski and De Raedq . l2003|) . Of course, 
not only the fermionic density matrix F but also its 
interacting counterpart, exp(— (3H), can be expanded 
in polynomials, which leads to similar methods for 
intera cting quantum systems ijlitaka and Ebisuzaki 
2003). 

To give an impression, in Figure 1171 we compare the 
optical conductivity of the Anderson model calculated 
with K PM (see Sec . IIII.D.2|) and with a projection ap- 
proach l|Iitakal . Il998[) . Over a wide frequency range the 
data agrees very well, but at low frequency the projec- 
tion results deviate from both KPM and the analytically 
expected power law cr(u>) — <jq ~ uja - Presumably this 
discrepancy is due to an insufficient resolution or a too 
short time-integration interval. There is no fundamental 
reason for the projection approach to fail here. 

In summary, the projection methods have a similarly 
broad application range as KPM, and can also compete 
in terms of numerical effort and computation time. For 
finite-temperature dynamical correlations the projection 
methods are characterized by a smaller memory con- 
sumption. However, in contrast to KPM they require a 
new simulation for each change in temperature or chemi- 
cal potential, which represents their major disadvantage. 

VI. CONCLUSIONS & OUTLOOK 

In this review we gave a detailed introduction to the 
Kernel Polynomial Method, a numerical approach that 
on the basis of Chebyshev expansion allows for an effi- 
cient calculation of the spectral properties of large matri- 
ces and of the static and dynamic correlation functions, 
which depend on them. The method has a wide range 
of applications in different areas of physics and quan- 
tum chemistry, and we illustrated its capability with nu- 
merous examples from solid state physics, which covered 
such diverse topics as non-interacting electrons in dis- 
ordered media, quantum spin models, or strongly cor- 
related electron-phonon systems. Many of the consid- 
ered quantities are hardly accessible with other meth- 
ods, or could previously be studied only on smaller sys- 
tems. Comparing with alternative numerical approaches, 
we demonstrated the advantages of KPM measured in 
terms of general applicability, speed, resource consump- 
tion, algorithmic simplicity and accuracy of the results. 

Apart from further direct applications of the KPM out- 
side the fields of solid state physics and quantum chem- 
istry, we think that the combination of KPM with other 



numerical techniques will become one of the major future 
research directions. Certainly not only classical MC sim- 
ulations a nd CPT, but poten tially also other cluster ap- 
proaches l)Maier et all 12005^ or quantum MC can profit 
from the concepts outlined in this review. 
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